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ABSTRACT 

Research  has  been  carried  out  on  observations  of  surface 
and  body  waves  from  explosions  and  earthquakes,  array  proces¬ 
sing  of  seismic  data,  source  properties  of  earthquakes,  theore¬ 
tical  and  experimental  studies  of  body  wave  propagation,  and 
structure  and  dynamics  of  the  crust  and  upper  mantle. 

The  spectra  of  waves  generated  by  earthquakes  and  under¬ 
ground  explosions  have  been  studied  using  the  Large  Aperture 
Seismic  Array  (LASA)  in  Montana  and  the  M.I.T.  tiltmeter  in 
Harvard,  Massachusetts. 

Shear  waves  generated  by  large  underground  explosions 
have  been  used  to  estimate  the  magnitude  of  tectonic  strain 
release  accompanying  the  explosions. 

Near-field  observations  of  the  Parkfield  earthquakes  are 

fitted  by  a  gliding-edge  dislocation  model.  Seismic  data  were 

used  to  investigate  the  focal  depth  of  earthquakes  and  core 
waves . 

Theoretical  studies  of  body  wave  propagation  have  been 
done  using  finite  differences,  complex  frequency,  and  genera¬ 
lized  ray  theory.  Non-linear  loss  mechanisms  are  being  treated 
numerically. 

The  structure  of  the  crust  and  upper  mantle  has  been 
studied  using  seismic  travel-times  and  spectral  ratios.  Nu¬ 
merical  simulations  have  been  done  on  the  thermal  behavior  of 


a  downgoing  slab  and  on  the  thermal  and  viscous  flow  benavior 
at  mid-ocean  ridges. 
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I. 


INTRODUCTION 


This  report  summarizes  research  carried  out  under  the 
Post-Doctoral  Program  in  Seismology  during  the  period  1  July 
1969  to  30  June  1971.  Staff  members  of  the  Department  of 
Earth  and  Planetary  Sciences  and  Lincoln  Laboratory  partici¬ 
pated  in  this  program. 

The  work  described  here  is  divided  into  five  categories: 

1)  observations  of  surface  and  body  waves  from  earthquakes 
and  explosions,  2)  properties  of  earthquake  sources,  3)  array 
processing  of  seismic  data,  4)  theoretical  and  experimental 
studies  of  body  wave  propagation,  and  5)  studies  of  the  struc¬ 
ture  and  dynamical  behavior  of  the  crust  and  upper  mantle. 

Rayleigh  waves  from  underground  nuclear  explosions  and 
earthquakes  in  western  North  America  recorded  by  a  mercury 
tube  tiltmeter  and  a  conventional  long-period  vertical  seismo¬ 
meter  at  Harvard,  Massachusetts,  have  been  used  to  investigate 
possible  spectral  ratio  techniques  for  discrimination  between 
explosions  and  earthquakes.  The  ratio  of  the  spectral  ampli¬ 
tudes  in  the  bands  from  15  to  22  sec  and  22  to  60  sec  is 
found  to  provide  complete  discrimination  for  events  in  the 
western  United  States,  but  not  for  the  Milrow  explosion  in  the 
Aleutians.  Regional  source  effects  appear  to  significantly  in¬ 
fluence  spectral  ratios  and  should  be  studied  further. 

Frequency-wave  number  spectra  of  Rayleigh  waves  from  ex- 


plosions  and  earthquakes  in  the  western  United  States  have 
been  studied  using  the  Montana  LASA.  A  loss  of  coherent  sig¬ 
nals  for  periods  near  50  sec  has  been  found  for  the  NTS 
explosions,  but  not  for  earthquakes  or  the  Colorado  explosion, 
Rulison.  This  phenomenon  is  not  well  understood,  but  may  be 
related  to  the  shallow  source  depth  of  NTS  explosions.  Various 
tests  have  indicated  that  these  observations  are  probably  not 
significantly  affected  by  non-linear  behavior  of  the  instru¬ 
ments  for  large  signals. 

Frequency-wave  number  spectra  of  Rayleigh  waves  from 
the  Milrow  explosion  show  that  significant  'multipathing' 
occurred  at  periods  in  the  range  20  to  25  sec.  No  'multipathing' 
was  observed  at  longer  periods . 

Shear  waves  generated  by  the  large  underground  nuclear 
explosions  Greeley,  Boxcar,  and  Benham  have  been  studied  with 
the  aim  of  determining  the  characteristics  of  strain  release 
in  the  focal  region.  Shear  waves  alone  do  not  furnish  suffi¬ 
cient  information  to  uniquely  specify  the  strength  and  orienta¬ 
tion  of  a  double-couple  source  assumed  to  coincide  with  the 
explosions;  two  of  four  parameters  must  be  determined  from  other 
types  of  observations.  If  this  is  done  by  assuming  that  fault¬ 
ing  is  of  the  vertical  strike-slip  type,  it  is  found  that  the 
strength  of  the  double-couple,  relative  to  the  explosion,  is 

0.75  +  0.06  for  Greeley,  0.20  +  0.02  for  Boxcar,  and  0.44  + 

0.04  for  Benham. 
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For  an  interpretation  of  near-field  seismic  motions,  a 
propagating,  vertical  strike-slip  fault  may  be  modeled  approx¬ 
imately  as  a  uniformly-gliding  edge  dislocation.  In  such  a 
model,  the  perpendicular  and  parallel  components  of  motion  at 
the  fault's  surface  are  simply  related.  Further,  the  amplitude 
of  the  displacement  is  sensitive  to  the  ratios  of  rupture-to- 
shear  wave  velocity  and  shear-to-compressional  wave  velocity. 

A  layer  of  low-velocity  sediments,  for  instance,  may  produce 
a  large  effect.  The  fault  model  is  applied  successfully  to 
near-field  observations  of  the  1966  Parkfield,  California, 
earthquake . 

The  focal  depth  of  the  three  largest  earthquakes  in  the 
Parkfield,  California,  sequence  was  determined  using  first  and 
second  P  arrivals  at  Berkeley,  California.  Data  were  insuf¬ 
ficient  to  determine  whether  the  second  arrival  corresponds  to 
refraction  along  a  discontinuity  within  the  crust  or  to  a  com¬ 
plicated  source-time  function.  With  either  interpretation, 
however,  the  focal  depths  for  the  three  earthquakes  apparently 
increased  with  time  through  the  sequence . 

A  study  of  core  waves  shows  that  reflection  at  a  discon¬ 
tinuity  in  the  outer  core,  near  the  inner-core  boundary,  pro¬ 
duces  precursors  to  PKIKP  observed  between  130°  and  142°. 
Neither  reflections  in  the  outer-core  nor  higher-mode,  dif¬ 
fracted  waves  can  produce  the  tail  of  the  diffracted  P  waves 


observed  between  105’  and  125“ .  It  appears  that  reflections 
or  multiple  paths  in  the  upper  mantle  must  be  responsible. 

Two  independent  methods  for  solving  wave-progation  pro¬ 
blems  in  heterogeneous  media,  numerical  integration  of  equa¬ 
tions  of  motion  using  finite-difference  techniques  and  complex- 
frequency  solutions  obtained  by  the  approximate  wave  theoreti¬ 
cal  method  of  Aki  and  Larner  a 9 70) ,  are  compared  for  several 
problems  and  are  found  to  be  in  excellent  agreement.  The  tech¬ 
niques  are  applied  to  the  response  of  soft  sedimentary  basins 
-o  shear  waves  and  underscore  the  need  to  account  realistically 
for  structural  heterogeneities  to  properly  estimate  surface 
motion  and  earthquake  risk  near  such  structures. 

The  propagation  of  acoustic  waves  in  a  fluid  medium  con¬ 
sisting  of  two  homogeneous  half-spaces  separated  by  a  linear 
transition  layer  has  been  studied  in  the  frequency  domain. 

The  solutions  obtained  are  exact  and  are  evaluated  by  numerical 
contour  integration  in  the  complex  wave  number  plane.  The 
solutions  are  compared  with  approximate  solutions  obtained  by 
the  saddle  point  method  and  the  sense  in  which  the  approximate 
solutions  are  asymptotic  ones  is  discussed.  The  contribution 

from  complex  poles  is  found  to  be  significant  for  waves  near 
the  critical  distance. 

The  extension  of  generalized  ray  theory  from  plane-parallel 
to  spherical  layers  is  straightforward  for  pulses  of  short 
duration  because  plane-wave  reflection  and  transmission  coeffi- 


cients  can  then  be  used.  Synthetic  seismograms  can  then  be 
calculated  using  the  Cagniard-de  Hoop  theory.  The  upper  limit 
of  pulse  duration  for  applicability  of  this  generalized  ray 
technique  to  body  wave  propagation  in  the  deep  mantle  is  75 
sec  for  S  waves  and  40  sec  for  P  waves.  The  technique 

has  been  used  to  obtain  preliminary  velocity  models  beneath 
North  America. 

Numerical  calculations  have  been  made  of  the  expected 
thermal  behavior  of  a  slab  of  lithosphere  as  it  sinks  into 
the  mantle,  as  is  thought  to  happen  beneath  island  arcs.  The 
relative  effects  of  lattice  and  radiative  thermal  conductivity 
and  viscous,  radioactive,  compressive  and  latent  heating  have 
been  investigated.  Seismic  wave  travel- times  and  amplitudes 
are  the  measurable  quantities  most  affected  by  slab  structure 
and  potentially  most  useful  for  further  refining  thermal  models. 
Waves  propagating  down  the  slab  from  shallow  earthquakes  are 
advanced  by  as  much  as  4  sec.  Further,  strong  focusing  and 
shadow-zone  effects  are  produced  in  the  radiation  pattern  of 
earthquakes  near  or  within  the  slab. 

Compressional  waves  observed  along  a  profile  extending 
from  the  western  United  States  to  the  Great  Lakes  have  been 
studied  to  resolve  details  of  the  structure  in  the  upper  1000 
km  of  the  mantle.  The  best  fitting  model  satisfies  the  ob¬ 
served  travel-times  and  waveforms.  Synthetic  seismograms 
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are  constructed  using  the  Cagniard-de  Hoop  method.  Large 
observed  later  arrivals  are  in  excellent  agreement  with  cal¬ 
culated  waves  related  to  discontinuities  near  430  km  and 
660  km  depths.  In  addition,  a  rather  high  velocity  gradient 
near  the  550  km  depth  is  needed  to  explain  large  arrivals 
at  about  21°  distance. 

The  crustal  structure  beneath  LASA  has  been  investigated 
using  the  spectral  ratio  of  vertical-to-horizontal  displace¬ 
ments  of  long-period  P  waves.  The  period  of  the  lowest  fre¬ 
quency  peak  in  the  spectral  ratio,  related  to  the  vertical  P 
wave  travel-time  in  the  crust,  is  quite  variable  over  the 
area  spanned  by  the  array.  Variations  in  crustal  thickness 
of  about  7  km  are  implied.  The  Moho  generally  shoals  from 
the  northeast  to  the  southwest  across  the  array  and  exhibits 
a  synclinal  structure  with  the  axis  plunging  toward  the  north¬ 
east  in  the  southwest  quadrant  of  LASA. 

A  numerical  simulation  of  upper  mantle  convection  has 
been  constructed  using  a  two-dimensional  time-dependent  model 
that  allows  large  viscosity  variations.  The  method  has  been 
applied  to  sea-floor  spreading  with  the  objective  of  examining 
the  driving  mechanism  of  mid-ocean  ridges.  Counterflow  below 
the  plates  is  confined  to  depths  less  than  340  km.  It  is 

| 

found  that  for  a  ^reading  velocity  of  1.2  cm/yr,  the  ridge 
can  produce  compressive  stress  in  the  lithosphere  out  to  a 
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distance  of  1600  km.  For  a  spreading  rate  of  6  cm/yr, 
however,  this  model  is  clearly  excluded  for  it  requries  an  ex¬ 
cessively  large  stress  in  the  lithosphere.  Therefore,  upwelling 
material  must  cross  the  seismic  discontinuity  at  the  400  km 
depth . 

Numerical  calcualtions  are  being  made  of  a  spherically 
symmetrical  wave  due  to  explosions  in  a  number  of  rock  types. 
Propagation  in  the  range  of  peak  stress  below  1  kb  is  cal¬ 
culated  starting  from  the  pulse  shape  calculations  of  earlier 
works.  An  anelastic  material  model  is  used  with  shear  hystere¬ 
sis  independent  of  the  rate,  corresponding  to  a  nominal  Q  for 
P  waves  of  180.  Preliminary  results  indicate  that  small 
irreversibilities  have  a  larger  effect  on  surface  waves  than 
on  body  waves. 

Some  of  the  work  discussed  herein  has  already  been  pub¬ 
lished.  The  summary  for  such  work  is  given  below  in  abstract 
form.  The  complete  results  may  be  found  by  consulting  the 
appropriate  reference  given  in  Section  VII  of  this  report. 


Reference 

Aki,  K.  and  K.  Larner,  Surface  motion  of  a  layered  medium  having 
an  irregular  interface  due  to  incident  plane  SH  waves, 

J.  Geophys .  Res.,  75,  933-954,  1970. 
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OBSERVATIONS  OF  SURFACE  AND  BODY  WAVES  FROM  EXPLOSIONS  AND 

EARTHQUAKES 

II. 1  Discriminations  of  Earthquakes  and  Explosions  by  the 

Rayleigh-Wave  Spectral  Ratio  by  John  S«  Derr  (Abstract) 

Limited  data  from  a  long- period  vertical  seismometer  and 
the  prototype  of  a  new  mercury  tiltmeter  at  the  station,  Har¬ 
vard,  Massachusetts,  tend  to  confirm  that  the  Rayleigh-wave 
spectral  ratio  of  a  short-  to  long-period  energy  provides  a 
useful  discriminant  for  underground  explosions  in  the  Western 
United  States  at  distances  of  from  28°  to  38°  in  the  magnitude 
range  of  3.5  <  Mg  <  5.6.  This  corresponds  to  discrimination 
in  the  range  of  5.0  <  mb  <  6.3.  The  same  spectral  ratio,  however, 
does  not  discriminate  the  large  Milrow  event  in  the  Aleutians 
from  earthquakes  in  the  same  region,  possibly  because  of  dif¬ 
ferences  in  earthquake  source  mechanisms  in  different  regions 
and  because  this  explosion  released  significant  tectonic  strain. 
The  Rayleigh  spectral  ratio  complements  the  Mg  -  m^  discriminant 
by  providing  a  clear  separation  of  earthquakes  and  explosions 
in  some  geographical  regions  where  earthquakes  are  of  low  stress 
drop  and  generally  strike-slip  in  mechanism. 

II.  2  Radiation  Patterns  of  S  Waves  from  Underground  Nuclear 
Explosions  by  Tomowo  Hirasawa  (Abstract) 

Use  is  made  of  polarization  angles  of  S  waves  from  the 


Greeley ,  Boxcar ,  and  Benham  events  to  determine  focal  mechanisms 
of  earthquakes  presumably  triggered  almost  instantaneously  by 
the  events.  Although  other  possibilities  are  r^.t  completely 
excluded ,  the  observed  S  waves  are  interpreted  as  mixed  waves 
direct  S  phase  and  reflected  phases  (sS,pS)  from  a  composite 
source  of  explosion  and  double-couple  force.  Assuming  a  verti¬ 
cal  strike  slip,  the  strike  of  the  P-wave  nodal  plane  and  the 
relative  strength  of  the  double-couple  source  compared  with 
that  of  the  explosion  are  determined.  The  seismic  moments  of 
triggered  earthquakes  are  calculated  as  3.9  ~  6.3  x  1023  dyne 
cm  for  Boxcar  and  1.4  ^  2.3  x  1024  dyne  cm  for  Benham.  Par¬ 
ticularly  in  the  case  of  Benham,  the  rise  time  of  about  2  sec 
together  with  the  seismic  moment  of  2.3  x  1024  dyne  cm  is  com¬ 
patible  with  the  observations  of  seismic  waves.  The  possibility 
of  estimating  the  pre-existing  tectonic  shear  stress  is  suggest¬ 
ed  in  the  case  of  Boxcar.  The  tectonic  stress  is  estimated  at 
130  bars  from  the  assumption  of  stress  relaxation  due  to  a 
cavity  created  by  the  nuclear  detonation. 


ARRAY  PROCESSING  OF  SEISMIC  DATA 


m .  1  Twenty-  to  Eighty-Second  Period  Characteristics  of 
Nuclear  Explosions  Recorded  at  LASA  by  Harry  Mack 

Abstract 

Surface  waves  recorded  at  LASA  generated  by  several  NTS 
events  have  been  subjected  to  frequency-wavenumber  analysis. 

All  the  spectra  show  a  loss  of  signal  at  0.02  Hz  but  there 
appears  to  exist  coherent  propagation  at  lower  frequencies, 
producing  a  notch  in  the  normalized  peak  wavenumber  spectrum 
at  0.02  Hz.  Several  earthquakes  and  the  deeply-buried  Colorado 
nuclear  explosion  RULISON  do  not  exhibit  this  notch.  This  may 
indicate  that  the  frequency-wavenumber  characteristics  peculiar 
to  NTS  events  is  a  result  of  shallow  source  depth. 


Frequency-wavenumber  (F-K)  spectra  have  been  calculated 
for  the  vertical  components  of  several  NTS  events  recorded  at 
LASA.  For  explosions  in  the  low-intermediate  yield  range  the 
wavenumber  spectra  exhibit  common  features.  There  is  a  loss  of 
coherent  signal  at  0.02  Hz  but  there  is  apparent  coherent  pro¬ 
pagation  at  lower  frequencies. 

Figures  la,  lb,  and  lc  show  high-resolution  F-K  spectra 
(Capon,  1969)  for  the  event  LANPHER  at  frequencies  of  0.05  Hz. 
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0.02  Hz,  and  0.015  Hz.  The  loss  of  signal  at  0.02  Hz  is  ob¬ 
vious  from  the  lack  of  a  coherent  peak  in  the  wavenumber  plane 
(Figure  lb) .  The  return  of  coherent  propagation  at  an  appar¬ 
ently  lower  frequency  is  illustrated  by  the  peak  in  Figure  lc 
which  has  a  phase  velocity  of  4.2  km/sec  which  is  about  right 
for  Rayleigh  waves  at  this  period.  The  apparent  direction  of 
propagation  is  further  west  than  the  true  back  azimuth.  However, 
the  resolution  of  the  array  is  poor  at  these  low  frequencies 
and  there  may  be  quite  a  large  error  in  velocity  shown  in  Fi¬ 
gure  2  for  LANPHER.  The  average  power  is  the  power  spectrum 
averaged  over  the  individual  sensors  of  LASA  and  has  contri¬ 
butions  from  both  coherent  and  incoherent  disturbances.  The 
values  are  normalized  to  the  power  at  a  frequency  of  0.05  Hz. 

There  is  an  increase  in  power  at  frequencies  lower  than 
0.02  Hz.  The  lower  curve  in  the  same  figure,  associated  with 
the  right-hand  scale,  is  the  normalized  peak  value  of  the  wave- 
number  spectrum  at  each  frequency  and  is  a  measure  of  the  ratio 
of  coherent  power  to  total  power  as  determined  by  F-K  analysis. 

A  distinct  drop  is  seen  in  this  value  at  frequencies  around 
0.02  Hz  corresponding  to  the  absence  of  coherent  propagation 
illustrated  in  Figure  lb.  The  curve  then  rises  again  towards 
the  lower  frequencies  and  is  maximum  at  a  frequency  of  0.015 
Hz  agreeing  with  the  coherent  peak  in  Figure  lc. 

The  frequency-wavenumber  characteristics  of  LANPHER  have 
been  found  for  six  other  events  of  the  low-intermediate  yield 
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range.  The  body  wave  magnitude  range  (mb)  of  these  events 
was  5.2  to  5.8.  So  far,  this  apparent  energy  at  long  periods 
has  not  been  observed  at  LASA  for  smaller  events.  This,  how¬ 
ever,  may  be  a  signal-to-noise  problem.  For  illustration. 
Figures  3a  and  3b  show  the  average  power  and  peak  wavenumber 
spectra  for  the  events  SLED  and  ZAZA.  The  shapes  of  the  wave- 

number  spectra  are  seen  to  be  very  similar,  both  with  the  notch 
at  0.02  Hz. 

F-K  analysis  of  earthquakes  located  in  Nevada,  California, 
Baja,  California,  and  off  the  coast  of  Oregon  has  failed  to 
produce  a  wavenumber  spectrum  with  a  notch.  The  spread  of 
surface  wave  amplitudes  at  0.03  Hz  for  the  earthquakes  is  the 
same  as  for  the  NTS  events .  Figure  4  shows  the  average  power 
and  peak  wavenumber  spectrum  for  the  event  from  Baja,  California. 
Both  curves  are  smooth  compared  with  those  shown  in  the  pre¬ 
vious  figures.  The  wavenumber  spectrum  curves  for  the  earth¬ 
quakes  were  not  all  similar  but  none  exhibited  a  notch  at  any 
frequency. 

The  analysis  procedures  were  designed  such  that  frequency 
or  wavenumber  windowing  could  not  produce  a  null  at  0.02  Hz. 
Analysis  of  the  high  amplitude  calibration  recording  rules  out 
non- linearities  in  the  amplitude  range  observed.  There  is  still 
the  possibility  of  high  frequency  surface  waves  causing  spurious 
movements  of  the  long-period  seismometer  mass  (Berkhemer  and 
Schneider,  1964).  The  argument  is  that  these  surface  waves 
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would  propagate  across  the  array  and  produce  a  long-period 
transient  from  each  seismometer.  F-K  analysis  of  these  tran¬ 
sients  would  produce  apparent  coherent  propagation  at  low  fre¬ 
quencies  with  low  phase  velocity.  Some  measured  phase  veloci¬ 
ties  were  low  and  this  would  agree  with  the  above  arguments. 
However,  this  is  not  always  the  case  as  can  be  seen  from  the 
velocity  of  4.2  km/sec  shown  in  Figure  lc.  This  higher  veloci¬ 
ty  is  hard  to  reconcile  with  expected  phase  and  group  velocities 
of  short-period  surface  waves.  Also,  the  short-period  surface 
waves  from  the  earthquakes  produced  no  such  effect. 

The  nuclear  event  RULISON,  located  in  Colorado,  produced 
F-K  spectra  similar  to  the  earthquakes.  Figure  5  shows  the 
power  and  wavenumber  spectrum  for  RULISON.  Coherent  propagation 
exists  at  0.02  Hz  and  disappears  at  lower  frequencies.  Although 
RULISON  had  a  body  wave  magnitude  of  5.3,  the  surface  wave  dis¬ 
placement  at  LASA  was  greater  than  most  of  the  other  NTS  events 
because  the  detonation  site  was  nearer.  RULISON  was  buried 
more  deeply  than  is  normal  at  NTS.  This  fact,  coupled  with  the 
earthquake  spectra,  may  indicate  that  the  wavenumber  spectra 
peculiar  to  NTS  events,  whether  real  or  not,  is  a  result  of 
shallow  source  depth. 
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Figure  1.  High  resolution  F-K  spectra  for  the  event  LANPIIER 
recorded  at  LASA  at  frequencies  of  la)  0.05  Hz,  lb) 
0.02  Hz,  and  lc)  0.015  Hz.  The  contours  are  in  dB. 
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Figure  2.  Average  power  spectrum  of  LANPIIER  recorded  at  I.ASA 
(upper  curve  and  left  hand  scale),  and  normalized 
peak  values  of  wavenumber  spectra  versus  frequency 
(lower  curve  and  right  hand  scale). 
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Figure  3a, 


Average  power  spectra  of  SLI£D  recorded  at  LASA 
(upper  curve  and  left  hand  scale),  and  normalized 
peak  values  of  wavenumber  spectra  versus  frequency 
(lower  curve  and  right  hand  scale). 
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Figure  5.  Average  power  spectrum  of  RULISON  recorded  at  LASA 
(upper  curve  and  left  hand  scale),  and  normalized 
peak  values  of  wavenumber  spectra  versus  f requcncy 
(lower  curve  and  right  hand  scale). 
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III. 2  Multipa thing  of  Rayleigh  Waves  Generated  bv  Milrow 
by  Harry  Mack 

Abstract 

Frequency  wavenumber  analysis  of  Rayleigh  waves  gneerated 
by  MILROW  shows  that  multipathing  occurs  between  the  source 
and  LASA  in  the  period  range  20-25  seconds.  No  multipathing 
is  detected  at  longer  periods. 


The  underground  nuclear  explosion  MILROW  was  detonated  on 
October  2,  1969  in  the  Aleutian  Islands.  The  body  wave  magni¬ 
tude  (n^)  was  6.5  and  Rayleigh  waves  were  well  recorded  at 
LASA.  Frequency-wavenumber  analysis  of  the  Rayleigh  wavetrain 
showed  that  significant  mutipathing  occurred  at  periods  in  the 
range  20-25  seconds  but  none  was  observed  at  longer  periods. 

Figure  1  is  a  composite  of  three  plots.  The  upper  curve 
shows  the  peak  wavenumber  spectral  value  versus  frequency. 

The  coherency  remains  high  and  almost  flat  from  0.06  Hz  to 
about  0.02  Hz.  The  values  then  decrease  rapidly  towards  the 
lower  frequencies.  The  triangle  shows  the  peak  wavenumber  spec 
tral  value  at  0.05  Hz  for  an  arrival  which  came  approximately 
five  minutes  after  the  original  wavetrain. 

The  middle  curve  shows  how  the  back  azimuth  of  arrival 
varied  with  frequency.  The  angle  increases  with  decreasing 


period  to  an  asymptotic  value  of  306  •  which  is  close  to  the 
true  back  azimuth.  The  triangle  indicates  that  the  later 
arrival  had  back  azimuth  of  321  •.  The  value  of  300-  at  0.06 
Hr  is  difficult  to  explain  except  that  this  value  may  be  the 
result  of  averaging  over  more  than  one  direction  of  arrival. 

The  bottom  curve  shows  the  phase  velocity  dispersion  curve. 
The  velocity  increases  with  period  in  an  expected  manner  andd 
then  low  values  appear  as  the  signal  disappears  at  the  lower 
frequencies .  This  phenomenon  is  the  result  of  spectral  leakage. 
Figure  2  shows  the  wavenumber  spectrum  at  0.05  Hz  at  the 

two  time  periods  indicated  on  the  figure.  The  main  arrival  at 
this  frequency  is  followed  by  ,  later  arrlval  ^  ^ 

wards  at  a  higher  angle.  Both  phase  velocities  are  the  same 

and  this  figure  illustrates  a  case  of  multipathing. 

in  conclusion,  the  above  results  indicate  that  multipathing 

of  Rayleigh  waves  may  not  be  a  problem  at  periods  greater  than 
40  seconds. 
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Figure  1.  Analysis  of  Rayleigh  waves  generated  by  MILROW  recorded 
at  LASA.  The  top  curve  is  the  normalized  peak  value  of 
the  wavenumber  spectra  versus  frequency.  The  middle  curve 
shows  direction  of  arrival  versus  frequency.  The  bottom 
is  the  phase  velocity  dispersion  curve.  The  triangles 
relate  to  the  later  arrival. 
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IV. 

SOURCE  PROPERTIES  OF  EARTHQUAKES 

IV. 1  A  Two-Dimensional  Moving  Dislocation  Model  for  a 

t  Abstract}**  Fault  by  D*  M‘  Z00™'  Aki,  and  T.  Todd 

For  a  propagating,  vertical,  strike-slip  fault  whose  break¬ 
age  extends  to  the  Earth's  surface,  previous  studies  by  Aki 
(1968)  and  Haskell  (1969)  have  suggested  that  the  near-field 
motions  may  be  similar  to  those  from  uniformly-gliding-edge 
dislocations.  The  theory  of  these  dislocations  in  uniform 
motion  leads  to  a  simple,  convenient  relation  between  the 
perpendicular  and  parallel  components  of  motion  at  the  fault's 
surface.  A  number  of  examples  are  considered  in  order  to  il¬ 
lustrate  this  relation  between  the  horizontal  components.  In 
general,  step-function-like  parallel  motions  result  in  pulse- 
like  perpendicular  displacements.  For  a  given  parallel  displace¬ 
ment,  the  amplitude  of  the  pulse  depends  in  a  concise  manner  on 
ratios  of  the  rupture  to  shear-wave  velocity,  and  on  shear-  to 
compressional-wave  velocity.  Increasing  either  of  these  ratios 
leads  to  an  increased  pulse-amplitude.  The  dislocation  model 
is  applied  to  the  near-field  observations  of  the  Parkfield 
earthquake.  The  resulting  estimate  of  the  total  fault  offset 
is  within  the  range  of  those  based  on  the  more  detailed  models 
of  Aki  (1968)  and  Haskell  (1969). 
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IV. 2  Focal  Depths  of  the  1966  Parkfield,  California, 
Earthquake  by  William  H.  Bakun 

Abstract 

Differences  in  arrival  times  of  seismic  phases  at  Berkeley 
California  (BRK) ,  A  =  270  km,  for  the  0408  GMT  June  28  (M  = 
5.1),  0426  GMT  June  28  (M  =  5.5),  and  1953  GMT  June  29  (M  =  5.0) 
1966  Parkfield,  California,  earthquakes  imply  that  focal  depths 
for  these  three  largest  events  of  the  1966  Parkfield  sequence 
increased  with  time  through  the  sequence.  The  data  available 
are  not  sufficient  to  determine  whether  the  observed  secondary 
arrivals  at  BRK  result  from  a  slower  propagation  path  or  are 
part  of  a  complicated  source- time  function. 


The  1966  Parkfield,  California  earthquake  sequence  occurred  on  the 
San  .Andreas  fault  in  central  California*  approicinatcly  equidistant  fro— 

San  Francisco  and  Los  Angeles.  The-  largest  earthquake  in  tne  sequence, 
at  01*26  GZ'IT  on  June  28,  1966,  had  a  Wood- Anderson  nagnitude  of  5.5. 

Surface  fracture  (Broun  and  Voider,  1967)  along  a  38  In  segment  o_  the 
San  Andreas  fault  zone  in  the  epicentral  region  accompanied  the  sequence. 

An  inport  ant  aspect  of  the  Parkfield  earthquakes  is  the  wealth  of  near¬ 
field  seismic  data  available  to  model  the  earthquake  source  and  uo  test 
techniques  of  estimating  earthquake  source  parameters  from  far-field 
seismic  data.  Consequently  the  1966  Parkfield  sequence  is  unique  —  at 
least  among  magnitude  5.5  sequences  --in  the  quantity  01  research  thau 
has  been  reported  concerning  its  focal  mechanism  and  source  parameters 
(AkL,  1967;  Eaton,  1967;  Filson  and  2*cEvilly,  1967;  foEvilly  gtjcl. , 

1967;  Aki,  1968;  Smith  and  Ifyss,  1968;  Ifa,  1968;  liyss  and  Brunc,  1968; 
Haskell,  1969;  Scholz  et  al. ,  1969;  Bcora  et  al.,  1970;  Eaton  et  al», 

1970). 

Eaton  et  al.,  (1970)  reported  that  9$%  of  the  Parkfield  aftershocks 

considered  in  their  study  had  focal  depths  between  1  and  12  km,  and  none 

were  deeper  than  15  km.  Available  focal  depths  for  large  aftershocks 

(M  ?  2.0)  in  the  sequence  range  from  1.6  to  12.2  km  with  an  average  depth 

near  5  kn  (1-IcEvilly  et  al..  1967).  Although  fecal  depths  are  available 

for  some  of  the  aftershocks,  near  data  are  not  sufficient  to  assign 

focal  depths  to  the  3  largest  events  (Table  1)  in  the  sequence.  The 

USCGS  lists  focal  depths  of  5,  U,  and  5  to®  for  the 

OI1O8,  01*26,  and  1953  C-I'xT  Parkfield  events  respectively  (USCGS  Earthquake 

Data  Report  PIE  #1*5-66).  The  relatively  large  standard  errors  assigned  to 

« 

the  focal  depths  (-6  tea)  by  the  USCGS  indicate  significant  scatter  ir.  the 
data.  Using  Love  and  Rayleigh  wave  spectral  amplitudes. 
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Tsai  and  AkL  (1970 )  estimated  *  depths  of  7,  k  csd  9  tea  respectively 
for  those  3  events.  In  this  report  seismograms  written  at  the  University 
of  California  seismographic  station  BRK  (A  .  2?0km)  are  interpreted  to 
imply  that  focal  depths  for  these  3  Park  field  events  increased  with  time 
through  the  sequence.  (It  should  be  stressed  that  this  conclusion 
regarding  focal  depth  is  applicable  only  to  these  3  largest  Parkfield 
events  and  not  to  the  sequence  as  a  whole.)  It  is  hoped  that  those  data 

will  serve  as  a  useful  constraint  to  discussions  of  the  source  parameters 
of  the  Parkfield  earthquakes. 


The  seismographic  station  BEK  (37°52.2*»  N,  122°l5.6«  tt,  elev.  -  8l  n) 
is  located  in  the  Havilland  Hall  vault  on  the  Berkeley  campus  of  the 
University  of  California.  The  broad  hand  3-component  system  is  described 
in  detail  by  Hodgers  etal.,  (1955).  Seismograms  written  at  BRK  for  the 
3  largest  1966  Parkfield  earthqualces  are  sho:n  in  figure  1.  The  traces 
are  aligned  at  the  trough  of  the  first  arrival,  labeled  Pn.  The  second 
arrival,  labeled  P2,  is  clearly  visible  on  each  of  the  traces.  The 
central  point  to  be  taken  from  figure  1  is  that  the  P2  -  Pn  time  ciffer- 
ence  increases  through  the  sequence  for  these  3  events.  The  absolute 
P1  ’  Pn  tine  ^^erence  depends  on  whether  the  time  picks  are  lead  from 
the  breaks  or  peaks  or  troughs  of  the  traces.  The  ?2  -  P^  tine  differ¬ 
ences  are  estimated  to  bo  2.9  sec  (0l*C8  fflff),  3.3  sec  (01*26  CBS)  and 
3.5  sec  (1953  GCT)  with  an  uncertainty  of  about  0.1  -  0.2  sec. 

The  3  Parkfield  events  shown  have  essentially  identical  epicenters  — 
the  difference  in  the  epicentral  coordinates  listed  in  Table  1  is  due 
to  the  reading  accuracy  <+  0.1  sec)  used  in  their  location.  The  epi¬ 
central  distance  to  BRK  is  about  270  km.  Eaton  (1967)  determined  an 


average  rupture  velocity  of  2.2  kn/sec  to  the  southeast  along  the  San 
Andreas  fault  for  the  01*2 6  GI2  event.  lilson  and  IfcEvilly  (1967)  found 
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a  rapture  velocity  of  2.2  ta/sec  and  a  fault  length  of  30  km  to  be  consis- 

tont  with  love  wave  spectral  amplitudes  at  BEX  for  the  0U2o  G®  event. 

Similarity  in  the  shape  of  the  love  wave  spectral  amplitudes  at  BEX  for 

tho  0^08  GIT  and  1953  GIT  events  and  2  smaller  1955  Par-field  earthquakes 

suggested  source  dimensions  loss  than  10  km  for  these  events  (?ilson  and 
IIcEvilly,  1967). 

Tho  fact  that  the  ?1  phaso  in  Figure  1  is  prominent  for  each  of  these 
3  apparently  quite  different  Paricfield  earthquake  sources  indicates  1) 
that  the  P1  phase  is  a  secondary  crastal  arrival,  independent  of  the  source, 
time  function  or  2)  that  ?1  is  a  source  phase,  common  to  the  3  Paricfield 
earthquakes,  possibly  the  'breakout  phase"  predicted  by  Savage  (1955)  end 
Bunldge  and  Ealliday  (1971).  It  is  unfortunate  that  the  seismic  data 
necessary  to  resolve  this  question  are  not  available.  Seismograph  sta¬ 
tions  at  near  and  rsgional  distances  in  California  were  generally  saturated 
by  the  first  P-wava  arrivals  from  these  3  events.  Records  from  Ilbod- 
Anderson  seismographs  operated  by  the  University  of  California,  Berkeley, 
at  Berkeley  and  1-bunt  Hamilton,  and  by  the  California  Institute  of  Tech¬ 
nology  at  Cottonwood,  Pasadena,  S^ta  Barbara  and  Tinemaha  are  riot 
saturated  but  the  horizontal  hbed- Person  instruments  are  not  suited 
for  recording  crustal  P-wave  arrivals.  (?or  sample,  tho  ?±  phase  re- 
coiMcd  on  the  broad-band  vertical  instrument  at  Berkeley  shown  in  Figure 
1  is  not  evident  on  the  Ivbod-Anderson  records  written  at  Berkeley).  At 
teleseisnic  distances,  the  ?-waves  from  these  Paricfield  events  are 
emergent  so  that  little  information  ccnceivmmg  secondary  arrivals  in 
the  P-wave  train  can  be  obtained  from  seismograms  written  at  telcseisr.de 
distances.  It  should  be  noted  that  the  standard  deviations  of  ?-::ave 
arrival  times  for  these  events  reported  in*  the  Bulletin  of  the  I-erna- 


tional  Seicnological  Centre  are  unusually  large  for  magnitude  $  earth¬ 
quakes#  One  explanation  of  those  large  deviations  is  that  stations 
reported  P-wave  readings  frees  different  phases  of  complicated  sources. 
Such  an  explanation  would  support  the  source  phase  hypothesis  for  ?^. 

In  this  report  the  P^  phase  is  considered  first  to  be  a  secondary  crustal 
arrival,  and  second  as  a  source  phase.  3h  each  case  ih3  result  that  the 
focal  depths  for  these  3  largest  Parlcfield  events  increased  with  time 
through  the  sequence  is  obtained. 
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P1  "  A  Secondary  Crustal  -Arrival 

o  interpret  the  P2  -  Pn  tame  differences  in  tones  of  different  focal 
depths  it  is  convenient  to  ess**  ray  paths  through  the  crust  and  toper 
mantle  for  the  Pn  and  ?1  phases.  P„  is  assumed  to  he  the  arrival  of 
a  head  nave  propagating  along  the  Kbho  and  ^  is  a  phase 

within  the  crust.  If  one  assumes  a  Pg  velocity  of  6  Wseo  and  the 

origin  tiir.es  given  in  i  +1.^  D  , , 

b  a*le  *hs  Pg  Ph£s®  (^rect  leave)  would  arrive 

at  BRK  about  8  sec  afte~  ■‘•he  P  .  , 

n  —rived.  A  sharp  arrival  near  this  tine 

is  clearly  evident  on  the  1*3  <K  trace  and  can  he  seen  for  the 
events  as  well.  *  Hgure  2  a  diagram  of  one  of  the  possible  ray  oaths 
of  the  h  phase  is  shown.  In  this  case,  Pj.  is  assumed  to  he  the  head 
wave  arrival  from  a  crustal  interface  at  a  depth  of  IS  to.  Assure 

the  crustal  model  and  ray  paths  shewn  in  Fi^ra  2,  .  Pn  tine  differences. 

computed  for  focal  depths  ranging  fro’p  2  to  K  ^  • 

°  ^  ro“’  2  t0  +5  km,  increase  with  focal 

depth.  Elis  result  —  that  p  d  +• 

%  the  P1  “  pn  tune  difference  increases  with 

focal  depth  -  is  not  affected  by  variations  in  the  assured  crustal  nodal. 

AS  en  enable  the  *  -  Pa  tine  differences  for  focal  depths  of  2  fa,  and 

W  bn  versus  the  oonpressional  velocity  in  the  bottom  layer  of  the 

crustal  node!  (6.8  Wseo  in  figure  2)  are  shown  in  figure  3.  »Ue  the 

cense  of  proportionality  between  the  P,  .  P„  tine  difference  and  focal 

depth  remains  the  sane,  it  is  clear  that  absolute  focal  depth  and  even 

absolute  focal  depth  difference  is  strongly  dependent  on  the  assured 

crustal  model*  Since  the  crustal  structure  ir.  central  California  is 

certainly  more  educated  than  the  model  assumed  in  these  calculations, 

precise  statements  concerning  the  focal  depths  for  the  3  largest  Pa-hiieM 

events  from  the  *  -  P„  data  prasanted  in  figure  1  are  not  justified. 


The  P-j,  - .  P  tine  differences  shown  in  Hgure  3  wore  computed  assuming 
that  the  ?n  and  phases  were  head  wave  arrivals  propagating  along  the 
I-bho  and  along  an  intermediate  crustal  interface  respectively.  Sinilar 
P-j^  *-  Pn  tine  difference  analyses  for  other  P  propagation  paths  within 
the  crust  lead  to  the  sane  result  —  that  an  increase  in  focal  depth 
results  in  an  increased  P^  -  Pn  time  difference.  In  fact,  it  is  clear 
that  the  P^  -  Pn  tine  difference  will  increase  with  focal  depth  as  long 
as  the  propagation  path  of  the  phase  is  confined,  to  the  crust  and 
the  Pn  phase  propagates  at  cr  near  the  top  of  the  mantle  at  a  greater- tha 
crustal  velocity.  The  conclusion  that  the  focal  depths  of  the  3  Parhfiel 
events  listed  in  Table  1  increased  with  time  in  the  sequence  then  is 
based  only  on  the  assumption  that  P^  is  a  crustal  phase  and  ?n  propagates 


at  or  near  the  top  of  the  mantle. 
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p  -  A  S01P.C3  3572CT 

Esrthqaalras  or.  the  San  -tadroas  fault  eppoe?  to  bo  well-suited  to 
tost  a  dynamic"  model  of  a  shallow  focus  earthquake  source  which  has  been 
proposed  by  toarUg,  and  Halllday  (1971).  &oir  'odsl  consists  of  a 
unifom  half-spaoo  that  is  prastrossed  to  proauoo  stroke- slop  r.or  .a* 
on  a  vertical  fault  plane.  Slip  is  initially  inhibited  on  the  nos- 
yielded  fault  surface  by  static  fsictior.  which  increases  with  depth. 

Initial  slip  results  in  a  stress  drop  which  drives  the  propagation  oe 
the  slip  upwards  and  downwards  iron  the  focus  at  the  shear  velocity/?. 

A  "breakout  phase"  is  generated  when  the  upwards  propagating  slip  reaches 

the  free  surface. 

Burridse  and  Halliday  (1971)  have  calculated  the  pulse  shapes  in 
the  far- field  radiation  and  the  residual  dit.xac^ents  and  stresses  on 
the  fault  plane  for  the  computationally  tractable  case  where  the  stress 
drop  on  the  fault  surface  decreases  quadratically  with  depth.  They  show 
that  in  this  case  displacement  on  the  fault  surface  extends  to  a  certain 
depth  which  is  independent  of  the  depth  of  focus.  This  surprises  con¬ 
clusion  is  in  accord  with  the  well-known  observation  that  earthquake  foca 
on  the  San  .Andreas  fault  do  not  extend  beyond  a  depth  of  about  IS  to. 
Burridse  and  Halliday  further  show  that  in  the  far-field  radiation,  for 
take-off  ansles  near  the  vertical,  the  "breakou-  ?--se  h-Cw.-.=  -w-° 
prominent  as  the  focal  depth  increases.  For  example,  if  the  calculations, 
of  Burridse  and  Halliday  are  scaled  so  that  the  depth  to  whica  dns?l-ce- 
asnt  on  the  fault  surface  extends  is  set  at  15  to,  then  the  "breakout 
phase"  is  larger  than  the  initial  motion  for  focal  depths  in  excess  of  b  to 
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A  possible  interpretation  of  the  se isnogroms  sho:r.  in  Figure  1  is 
that  the  P-^  phase  is  a  source  effect,  the  'breakout  phase"  predicted  by 

Bur  ridge  and  Halliday.  In  this  interpretation  P-^  propagates  from  the 

epicenter  to  the  Koho  as  a  P-uavc  and  then  as  a  head  rave  along  the 

Kbho.  The  -  Pn  tine  difference  is  thus  due  to  the  delay  of  ?- 

resulting  fro:n  propagation  from  the  fccus  to  the  free  surface  at  the 

shear  velocity  and  from  the  free  surface  doun  to  the  depth  of  focus  at 

the  comers  ssional  velocity.  It  is  clear  that  v.’ith  this  interpretation, 

an  increase  in  fecal  death  results  in  an  increased  Pi  -  P  time  differ- 

*  n 

ence  and  focal  depths  for  the  3  Par'.cfield  events  discussed  here,  icculd 
increase  from  about  7  km  for  the  OhOS  GIT  event  to  about  8  1cm  for  tho 
0li26  (HIT  event  to  about  9  1cm  for  the  1953  Ci-ZT  event. 
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DISCUSSIO’I 


<£v 


Tsai  and  Aid.  (1970)  determined  focal  depths  for  the  3  Parkfield  events 
studied  here  from  10-90  sec  love  and  Rayleigh  nave  spectral  amplitudes 
recorded  at  Alburquerque,  New  Nendo o  and  Colden,  Colorado.  Assuring  a 
point  source  vertical  strike-slip  dislocation,  they  found  focal  depths 
of  7  km  for  the  0*03  GIT  event  and  9  km  for  the  1953  C2  event  to  be 
consistent  with  the  observed  surface  wave  spectral  amplitudes.  Assuming 
a  vertical  strike-slip  dislocation  of  37  km  length  and  a  constant  rap¬ 
ture  velocity  of  2.2  kr/sec,  they  found  that  a  focal  depth  of  *  Jan  best 
fits  the  vertical  component  Rayleigh  wave  spectral  amplitude  recorded  at 
Atlanta,  Georgia  for  the  0*26  GKP  event.  Filson  and  NcSvilly  (l 967) 
and  \h  (1963)  have  suggested  that  the  0*26  GIT  event  could  have  been  a 
series  of  events  closedly  spaced  in  time.  The  variation  in  rupture 
velocity  resulting  from  a  multiple  source  would  complicate  the  surface 
wave  spectral  amplitudes  considerably  and  possibly  introduce  a  significant 
error  an  the  focal  depth  determined  for  the  0*26  GHT  event  by  Tsai  and 
ine  P]_  Pq  time  difference  cata  reported  here  support  the  focal 
depths  determined  by  Tsai  and  Aid.  in  that  their  estimates  of  focal  depths 
increase  in  time  through  the  sequence  for  the  0*08  GMT  and  1953  GIT 
events.  The  point  source  dislocation  model  assumed  for  these  2  events  is 
in  accord  with  the  love  wave  spectral  amplitude  data  of  PLlson  and 
NcSvilly  (1567).  Tsai  and  Aki»s  conclusion  that  the  focus  of  the  0*26 

GIIT  SVent  T,as  shallower  aan  the  "ooms  of  the  0*08  GMT  event  could  result 
from  errors  introduced  hy  ass  rains  too  simple  a  source  dislocation  model 

for  the  0*26  GIT  event.  If  the  dislocation  is  not  constant  but  varies 
with  depth,  then  Tsai  and  Aldis  results  are  more  indicative  of  an  average 
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location  depth  i.  n  the  depth  of  initial  rapture.  Regardless  of  the 
se  of  the  focal  depth  errors,  the  result  obtained  here  -  that  fecal 
ths  increased  -with  time  for  these  3  events  —  implies  an  empirical 
■or  of  about  for  focal  depths  less  than  10  kn  that  are  determined 
the  technique  proposed  by  Tsai  and  /id.  (1970). 
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TABLE  1 


1966  Parkfield,  California  earthquakes. 

Data  are  taken  from  1-IeEviIly  cl.  (1967 ) 


Date 


Origin  Time  (GMT)  Latitude  (N) 


Longitude  (\V) 


Magnitude  ( 


35°57. 6' 
35°57. 3* 
35°56. 6' 


120°30. 3' 
120°29.  9' 


5.1 


5. 5 


5. 0 


June  28,  1966 
June  28,  1966 
June  29,  1966 


04:08:56.2 

04:26:13.4 

19:53:25.9 


120°31. 5' 
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Figure  1 


Figure  2 


Figure  3 


FIGURE  CAPTIONS 


BRX  VERTICAL  seismograms  for  the  (top)  0h03  GMT, 
(middle)  01*26  GKP  and  (bottom)  1953  GI-TT  1956  Parkfield 


earthquakes.  The  traces  are  aligned  at  the  first  trough 


of  the  first  arrival,  labeled  ?n.  A  dashed  vertical 
line  has  been  drawn  through  the  first  trough  of  the  ?1 
phase  on  the  middle  seismogram  to  facilitate  compari¬ 
son.  Compression  is  devm  on  the  records. 


Diagram  (not  to  scale )  of  the  ray  paths  assumed  for  the 
Pn  and  P-l  phases  for  the  computation  of  the  P-^  -  pn 
time  differences  plotted  in  Figure  3.  The  crustal  model 
is  based  on  the  model  proposed  for  the  Parkfield  area  by 
Eaton  (1967). 


P1  “  pn  time  differences  for  focal  depths  of  2  km  and 

15  km  vs.  compressions!  velocity  at  the  base  of  the 

crust  (6.8  km/sec  in  Figure  2).  The  assumed  P  ’and  P, 

n  1 

ray  paths  and  crustal  model  shewn  in  Figure  2  were 
used  in  ohe  computations .  horizontal  arrows  indicate 
the  P}_  -  Pn  time  differences  taken  from  Figure  1. 
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V. 

THEORETICAL  AND  EXPERIMENTAL  STUDIES  OF  BODY-WAVE  PROPAGATION 

V.l  The  origin  of  precursors  to  Core  Waves  by  E.  Husebye 
and  R.  Madariaga  (Abstract) 

The  origin  of  the  precursors  of  the  core  waves  in  the 
range  105°  -  142°  is  studied.  Between  105°  and  125°  a  long 
tail  is  observed  after  the  P  wave  diffracted  by  the  core.  In 
the  range  130°  =  A  =  142°  we  usually  observe  short-period 
onsets  a  few  seconds  before  PKIKP;  these  are  the  waves  called 
P(GH).  Reflection  at  a  discontinuity  in  the  outer  core,  near 
the  inner-core  boundary,  is  shown  to  produce  the  P(GH)  branch. 
Reflections  in  the  outer  core  are  rejected  as  a  mechanism  for 
the  tail  of  the  P  diffracted  wave.  A  theoretical  study  of 
diffraction  of  P  by  the  core  shows  that  higher  modes  of  dif¬ 
fracted  waves  cannot  explain  the  observations  of  the  tail  of 
P  diffracted.  We  conclude,  by  elimination,  that  it  is  due  to 
reflections  or  multiple  paths  in  the  upper  mantle. 


Comparison  of  two  independent  Methods  for  the  Solution 
of  Wave  Scattering  Problems:  Response  of  a  Sedimentary 
Basin  to  Vertically  Incident  SH  Waves  by  D.M.  Boore, 
K.L.  Larner  and  K.  Aki  (Abstract) 


The  transient  solution  to  several  problems  that  was  obtained 
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by  numerical  integration  of  equations  of  motion  using  a 
finite  difference  (FD)  technique  is  compared  with  the  complex- 
frequency  solutions  obtained  by  the  approximate  wave  theo¬ 
retical  (AL)  method  of  K.  Aki  and  K.L.  Larner.  The  excellent 
agreement  between  the  two  solutions  not  only  provides  a  com- 
paritive  check  on  the  accuracies  of  the  two  techniques,  but 
also  demonstrates  that  the  interpretation  of  the  AL  solution 
is  comparable  to  the  Fourier  transform  of  the  transient 
solution  premultiplied  by  an  exponential  window.  Most  of  the 
paper  is  devoted  to  a  discussion  of  two  models  that  are  rele¬ 
vant  to  the  engineering-seismological  study  of  earthquake 
motions  in  soft  layers  of  varying  thicknesses.  The  FD  and  AL 
solutions  show  that  lateral  reverberations  of  waves  produced 
by  the  nonplanar  structure  form  complex  interference  patterns 
that  are  not  predicted  by  the  usual  flat-layer  approximations. 
In  one  example,  constructive  interference  enhances  the  peak 
amplitude  of  the  transient  motion  over  the  center  of  the  basin 
by  a  factor  of  3  relative  to  the  flat-layer  solutions.  The 
results  indicate  that  a  realistic  appraisal  of  earthquake 
hazards  in  areas  underlain  by  soft  surficial  layers  should 
include  the  effect  of  nonuniformity  in  the  structure. 
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V'3  ?r°">  3  Linear  Transition  Layer 

in  a  Fluid  Medium  by  T.  Hirasawa  and  M.J.  Berry  (Abstract) 

Reflected  and  head  waves  from  a  linear  transition  layer 

between  two  homogeneous  media  are  studied  in  the  frequency 

domain.  A  point  source  is  located  in  the  upr ;r  fluid  of  lower 

velocity,  and  the  velocity  structure  is  considered  to  be  con- 
tinuous  throughout. 

Exact  solutions  are  derived  by  numerical  evaluation  of 
the  contour  integrals  in  the  complex  wave-number  plane.  These 
are  compared  with  approximate  solutions  obtained  by  the  saddle 
point  method.  It  was  found  that  the  approximate  solutions  for 
head  and  reflected  waves  beyond  the  critical  distance  may  be 
regarded  as  the  asymptotic  ones  for  large  values  of  r/2h. 
Similarly,  the  approximate  solutions  for  reflected  waves  inside 
the  critical  distance  are  the  asymptotic  ones  for  large  values 
of  H/2h,  where  H  is  the  sum  of  vertical  distances  of  the  source 
and  a  receiver  from  the  transition  layer  of  thickness  2h,  and 
r  is  the  horizontal  distance  of  a  receiver  from  the  source. 

At  high  frequencies,  the  spectral  amplitudes  of  reflected 
waves  inside  the  critical  distance  are  proportional  to  w'1, 
while  head-wave  amplitudes  are  proportional  to  w"2/3  (u  belng 
the  angular  frequency) . 

Numerical  calculations  also  show  that  contributions  from 
complex  poles  are  significant  near  the  critical  distance  if 
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the  transition  layer  is  thick.  In  this  region,  the  amplitude 
variations  of  the  sum  of  the  reflected  and  secondary  waves 
are  similar  to  those  from  a  sharp  discontinuity,  although 
the  maximum  amplitude  for  a  wave  of  narrow-frequency  band 
width  occurs  at  greater  distances  for  thicker  transition  layers. 

V.4  Generalized  Ray  Theory  for  a  Layered  Sphere  by  F.  Gilbert* 
and  D.V.  Helmberger  (Abstract) 

Pulse  propagation  in  a  layered  sphere  can  be  investi¬ 
gated,  in  an  approximate  way,  by  what  has  come  to  be  known  as 
Cagniard's  method.  Classical  methods  are  used  in  the  analysis 
through  the  application  of  the  Debye  ray  expansion.  At  this 
stage  Lamb's  observation,  that  the  eikonel  is  linear  in  the 
frequency,  is  employed  to  bypass  the  usual  methods  for  evalu¬ 
ation  of  the  inverse  transform  integrals.  The  transient  re¬ 
sponse  for  each  Debye  ray  is  obtained  directly.  It  is  esti¬ 
mated  that  the  method  can  be  applied  to  mantle  S  pulses  with 
periods  less  than  75  s.,  and  to  mantle  P  pulses  with  periods 
less  than  40  s.  Preliminary  results  on  lateral  heterogeneity 
beneath  North  America  are  presented. 


*  At  the  Institute  of  Geophysics  and  Planetary  Physics,  Scripps 
Institution  of  Oceanography,  University  of  California  San 
Diego,  La  Jolla,  California  92037. 
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V.5  Calculated  Attenuation  of  an  underground  Explosion 

Wave  in  the  Nearly-Elastic  Range  by  D.J.  Andrews  and 
S.  Shlien. 


Abstract 

A  numerical  calculation  is  presented  of  a  spherically 
symmetrical  wave  in  granite.  Propagation  in  the  range  of 
peak  stress  below  1  kilobar  is  calculated,  starting  from 
the  pulse  shape  calculated  by  Butkovich  for  the  Hardhat 
event  at  100  meters  from  the  detonation  center.  An  anelastic 
material  model  is  used  with  shear  hysteresis,  independent 
of  the  rate,  corresponding  to  a  nominal  Q  for  P  waves  of  180. 
It  is  found  that  Fourier  components  of  different  frequencies 
do  not  attenuate  independently,  but  that  energy  is  transferred 
from  the  dominant  frequencies  around  15  Hertz  to  frequencies 
below  3  Hertz. 
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Observed  ground  shocks  from  underground  explosions 
in  the  close-in  region  have  pulse  lengths  much  shorter 
than  waves  observed  at  teleseismic  distances.  These  wave 
forms  have  been  matched  by  numerical  calculations  in  the 
highly  non-linear  and  irreversible  region  near  the  source 
(Butkovich,  1965) .  For  the  Hardhat  event,  5  kilotons  in 
granite,  Butkovich  calculated  a  peak  stress  at  100  meters 
of  1  kilobar  with  a  pulse  width  at  half  maximum  of  10 
millisec.  Beyond  this  radius  one  would  expect  the  granite 
medium  to  behave  nearly  elastically.  For  spherically  symmetric 
wave  propagation  in  the  elastic  region  a  component  would 
appear  with  wavelength  comparable  to  the  radius  of  the  ine¬ 
lastic  region  (Sharpe,  1938) .  No  transfer  of  erergy  to 
longer  wavelengths  would  occur  if  the  medium  behaved  elastically. 

Of  course  the  close-in  waveform  consists  of  Fourier 
components  of  all  frequencies.  Distant  waveforms  might  be 
explained  by  propagation  and  attenuation  of  each  Fourier 
component  independently  of  the  others,  as  in  the  case  of  a 
linear  viscoelastic  medium.  Then  the  energy  present  in  low 
frequencies  in  the  distant  waveform  is  a  remnant  of  the 
energy  present  in  those  components  of  the  close-in  waveform. 

The  mechanism  of  attenuation  in  crustal  rocks  has  been 
reviewed  by  Knopoff  (1964)  and  by  Gordon  and  Davis  (1968)  . 

They  find  that  in  polycrystalline  samples  the  stress-strain 


Preceding  page  blank 


52 


path  is  not  reversible/  but  has  a  small  hysteresis  that  is 
nearly  independent  of  strain  rate.  This  effect  is  not 
plasticity  in  the  usual  sense,  for  the  loading  and  un¬ 
loading  slopes  are  only  slightly  different.  The  effect 
probably  arises  from  friction  at  grain  boundaries. 

To  explain  the  change  of  waveform  from  close-in  to 
distant  stations  the  static  hysteresis  of  crustal  rocks 
should  be  taken  into  account.  For  this  type  of  irrevers¬ 
ibility  the  principle  of  superposition  does  not  hold,  so 
that  different  Fourier  components  will  not  propagate 
independently.  Since  propagation  in  such  a  medium  is  not 
amenable  to  analytic  solution,  a  numerical  finite  difference 
method  must  be  used.  A  method  of  Wilkins  (1964),  based  on 
direct  numerical  approximation  of  first  principle  equations, 
is  adopted. 

A  complete  prediction  calculation  of  the  distant  wave¬ 
form  would  have  to  account  for  the  free  surface  and  layers 
below  the  surface.  Such  a  calculation  would  have  to  be 
done  in  two  space  dimensions  with  axial  symmetry,  if  not 
in  three  dimensions,  and  would  require  rather  fine  zoning 
to  get  sufficient  accuracy.  Such  a  calculation  would  require 
a  large  amount  of  computer  time. 

In  this  work  a  spherically  symmetric  calculation  is 
done.  Although  it  is  not  a  complete  prediction  calculation. 
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it  shows  to  what  extent  the  assumption  of  superposition  of 
independent  components  is  in  error. 

The  material  chosen  for  the  calculation  is  reversible 
in  dilatation  and  irreversible  in  shear.  Shear  stress  is 
plotted  against  shear  strain  in  Figure  1  for  a  cycle  of 
loading  and  unloading.  The  initial  loading  slope  corresponds 
to  a  shear  modulus  of  0.306  megabar,  and  the  initial  un¬ 
loading  is  6  per  cent  steeper  with  a  shear  modulus  of 
0.326  megabar.  If  unloading  continues  until  shear  stress 
reverses  sign,  the  modulus  goes  back  to  the  smaller  value, 
and  then  again  becomes  steeper  for  reverse  unloading.  A 
closed  stress-strain  path  is  shown  in  Figure  1,  but  it  need 
not  be  closed,  depending  on  the  strain  history.  The  6  per 
cent  irreversibility  in  shear  corresponds  to  a  3  per  cent 
irreversibility  for  a  P  wave,  or  a  nominal  Q  of  180. 

Bulk  modulus  and  shear  modulus  chosen  for  the  calculation 
are  appropriate  to  granite  with  a  P  wave  velocity  of  5.4  km/sec. 
In  order  to  avoid  the  highly  irreversible  region  very  close  in, 
a  stress  pulse  was  applied  on  a  virtual  surface  at  a  radius  of 
100  meters.  Butkevich's  calculated  pulse  was  matched  roughly 
with  an  initial  radial  stress  of  1  kilobar  and  an  expo¬ 
nential  decay  with  a  decay  time  of  10  millisec.  If  stress 
decayed  to  a  finite  value  (due  to  gas  in  the  cavity)  rather 
than  to  zero,  low  frequency  components  would  be  larger. 


The 
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finite  difference  zone  size  used  in  the  calculation  is 
80  meters. 

A  pair  of  calculations  was  done;  one  with  a  reversible 
model  and  one  with  the  irreversible  model  shown  in  Figure  1. 
Waveforms  of  velocity  as  a  function  of  time  at  radii  of  0.9, 
1.7,  2.5  and  3.3  kilometers  were  computed  and  are  shown  at 
the  distant  stations  for  the  reversible  case  in  Figure  2, 
and  for  the  irreversible  case  in  Figure  3.  It  is  evident 
that  peak  velocity  is  smaller  in  the  irreversible  case. 

There  is  a  break  in  the  slope  of  the  waveform  in  the  irre¬ 
versible  case  when  the  stress  deviator  reverses  sign.  This 
is  due  to  the  sudden  change  of  shear  modulus  in  the  model 
and  should  not  be  considered  realistic.  An  important  fea¬ 
ture  not  evident  in  Figure  3  is  that  the  integral  of  the 
waveform  is  not  zero,  but  that  there  is  a  permanent  dis¬ 
placement.  The  permanent  displacement  depends  on  the 
amount  of  irreversibility, ‘but  not  on  details  of  the  hys¬ 
teresis  loop. 

Calculated  velocity  waveforms  at  the  4  stations  in 
each  calculation  were  Fourier  analyzed  by  the  Cooley-Tukey 
algorithm.  The  time  window  used  for  each  waveform  is  320 
milliseconds,  so  that  Fourier  amplitudes  are  found  at  fre¬ 
quency  intervals  of  3.12  hertz.  Velocity  spectra  are 
plotted  for  the  reversible  case  in  Figure  4  and  the  irre¬ 
versible  case  in  Figure  5.  Amplitudes  for  the  dominant 
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frequencies  of  15  to  20  hertz  are  smaller  for  the 
irreversible  case,  but  amplitudes  at  very  low  frequencies 
are  larger.  This  is  because  the  integrals  of  the  irre¬ 
versible  wave  forms  are  nonzero. 

In  seismological  literature  it  is  customary  to  measure 
attenuation  in  terms  of  the  coefficient  Q  given  by 


2  it  E 


where  E  and  AE  are  the  peak  energy  and  energy  loss  in 
a  single  cycle.  This  definition  implicitly  assumes  that  the 
medium  behaves  linearly,  such  that  the  ration  of  E  to  AE 
does  not  depend  upon  E  or  any  other  property  of  the  material 
which  is  affected  by  the  propagation  of  a  wave.  Therefore 
if  a  plane  wave  is  propagating  through  a  linear  medium,  then 
its  amplitude  should  decay  .exponentially  with  distance  given 


Uc  (  2c  Q  ) 


Q  may  be  a  function  of  frequency.  For  example,  for  a 
Kelvin  viscoelastic  material  Q  is  linearly  proportional 
.  In  the  case  of  earth,  both  laboratory 
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experiments  and  seismic  observations  show  that  Q  is 

independent  of  frequency.  Thus  high-frequency  seismic 

signals  decay  faster  than  low-frequency  signals  in  a  fixed 
distance. 


To  relate  the  finite  difference  solution  for  the  non- 

reversible  case  to  seismic  observations,  Q  „as  determined 

from  the  amplitude  spectrum  at  the  four  stations.  For  a 

spherically  symmetric  wave,  Blake  (1952)  showed  that  without 

any  attenuation  geometry  requires  a  harmonic  wave  to  propagate 
as 


i  /  \  /  .  /v  i  i  to  t  Ic  r  *  •'jO 

u  =  -j  Cu)A  f'C'1)  ^  I  4-  jU  <-/cJ  (L 


where  u  is  velocity,  r  radial  distance,  c  sound  speed, 

and  A'  and  p  constants.  Therefore  if  „e  include 
attenuation,  then 


U,  -  (  (  +  J^r/C  J  t  Vp  ( 


The  attenuation  factor,  1/a,  found  in  ;his  way  lg 
plotted  for  the  reversible  case  in  Figure  6.  The  attenu¬ 
ation  is  quite  close  to  that  expected  for  the  artificial 
viscosity  used.  This  is  a  check  on  the  accuracy  of  the 
numerical  method.  The  attenuation  from  artificial  vis¬ 
cosity  is  not  significant  at  frequencies  of  15  hertz  and 
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below.  If  the  principle  of  superposition  held,  one  would 
expect  that  1/Q,  would  be  increased  at  each  frequency  by 
1/180  or  0.0055  in  the  irreversible  case. 

The  attenuation  factor,  1/Q,  for  the  irreversible 
case  is  plotted  in  Figure  7  as  a  function  of  frequency  for 
three  pairs  of  stations.  The  three  curves  would  be  identical 
if  the  attenuation  were  exponential . 

Since  the  curves  definitely  do  not  coincide  it  implies 
that  Q  depends  upon  the  nature  of  the  signal  propagating 
through  the  medium.  For  instance,  the  Q  determined  would 
depend  on  whether  or  not  the  signal  is  dispersed.  There¬ 
fore  Q  is  no  longer  just  an  intrinsic  property  of  the 
medium.  Not  only  does  it  depend  on  the  frequency  for  which 
it  was  computed  but  also  On  the  presence  or  absence  of 
various  other  frequency  components  in  the  signal. 

Attenuation  at  the  dominant  frequencies  of  15-20  hertz 
is  larger  than  would  be  expected.  At  3  hertz  1/Q  is 
negative  and  is  off  the  scale  of  Figure  7.  Values  of  1/Q 
for  the  three  pairs  of  stations  are 


at 

-  .456 

0.9  -  1.7  km 

-.129 

1.7  -  2.5 

-  .022 

2.5  -  3.3 

at  3  hertz 
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Different  Fourier  components  are  not  propagatimg 
independently,  but  energy  is  being  transferred  from  the 
dominant  components  to  the  3  Hertz  component,  the  lowest 
frequency  component  that  can  be  analyzed  with  the  time 
window  used. 

In  the  work  that  is  in  progress  at  present,  the  reduced 
displacement  potential  is  being  examined  for  explosions  in  a 
number  of  rock  types.  In  an  elastic  material,  of  course,  the 
reduced  displacement  poteritial  is  invariant  as  the  wave  pro¬ 
pagates.  The  small  irreversibilities  used  in  numerical  cal¬ 
culations  have  a  significant  effect  on  the  reduced  displace¬ 
ment  potential.  Recall  that  displacement  in  surface  waves  is 
proportional  to  this  potential  rather  than  to  its  derivative 
as  is  the  case  for  body  whves.  Preliminary  results  indicated 
that  small  irreversibilities  have  a  larger  effect  on  surface 
waves  than  on  body  waves. 
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Figure  Captions 


FIGURE  1 


FIGURE  2 


FIGURE  3 


FIGURE  4 


FIGURE  5 


FIGURE  6 


FIGURE  7 


1  Shear  stress  versus  shear  strain  for  an  irre¬ 
versible  medium. 

2  Computed  radial  particle  velocity  as  a  function 
of  time  for  a  reversible  medium  at  distances  of 
2-5  and  3.3  kilometers  from  a  5  kton  explosion. 

I  Computed  radial  particle  velocity  as  a  function 
of  time  for  an  irreversible  medium  at  distances 
of  2.5  and  3.3  kilometers  from  a  kton  explosion. 
Square  root  of  energy  spectra  of  radial  particle 
velocity  for  a  reversible  medium  computed  at  four 
distances.  Units  of  spectra  are  in  cm/sec  -  millisec. 
Square  root  of  energy  spectra  of  radial  particle 
velocity  for  an  irreversible  medium  computed  at 

four  distances.  Units  of  spectra  are  in  cm/sec  - 
millisec. 

Q  versus  frequency  for  a  reversible  medium 

computed  at  3  pairs  of  adjacent  stations  using 
equation  [1] . 

Q  versus  frequency  for  an  irreversible  medium 
computed  at  3  pairs  of  adjacent  stations  using 
equation  [1J . 
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VI. 

STRUCTURE  AND  DYNAMICS  OF  THE  CRUST  AND  UPPER  MANTLE 

VI. 1  Temperature  Field  and  Geophysical  Effects  of  a  Down¬ 
going  Slab  by  M.N.  Toksflz ,  J.W.  Minear  and  B.R.  Julian 
(Abstract) 

The  factors  affecting  the  thermal  behavior  of  a  litho¬ 
spheric  slab  descending  into  the  mantle  are  so  numerous  and 
complicated  that  only  numerical  methods  can  accurately  account 
for  them.  A  series  of  finite-difference  calculations  of  the 
temperature  field,  including  the  effects  of  viscous  dissipation, 
adiabatic  compression,  phase  changes  and  radioactive  heat 
generation  are  carried  out,  and  the  relative  effects  of  dif¬ 
ferent  parameters  are  investigated.  An  analysis  of  the  stability 
and  convergence  of  the  numerical  method  indicates  that  the 
errors  are  small  and  can  be  reduced  to  any  desired  level  by 
varying  the  grid  size  and  the  time  steps.  At  a  crustal  spread¬ 
ing  rate  of  8  cm/yr,  with  all  heat  sources,  the  slab  reaches 
thermal  equilibrium  with  the  surrounding  mantle  at  a  depth  of 
about  650  km.  Among  observable  geophysical  quantities,  seismic 
travel  times  and  amplitudes  provide  the  most  information  about 
the  slab.  Surface  heat  flow  is  sensitive  to  subsurface  con¬ 
ditions  that  are  at  relatively  shallow  depths,  and  gravity 
anomalies  are  broad  and  are  masked  by  crustal  effects.  Three 
dimensional  calculations  predict  strong  bending  of  seismic 
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rays  near  slabs,  which  causes  strong  focusing  and  produces 
shadow  zones.  Analysis  of  travel-time  data  for  the  Tonga- 
Fiji  region  indicates  that  waves  propagating  down  the  slab 
from  shallow  events  are  advanced  by  about  4  sec.  Observa¬ 
tions  and  theoretical  travel-time  anomalies  based  on  calcu¬ 
lated  temperature  fields  are  in  general  agreement. 


VI. 2  Upper  Mantle  Structure  of  Midwestern  United  States  by 
D.  Helmberger  and  R.A.  Wiggins  (Abstract) 

Body-wave  observations  from  nuclear  events  at  the  Nevada 
Test  Site  and  several  earthquakes  near  the  western  edge  of 
the  United  States  have  been  used  to  construct  a  model  of  the 
upper  mantle  along  profiles  extending  toward  the  Great  Lakes. 
The  Cagniard-de  Hoop  technique  for  computing  synthetic  seis¬ 
mograms  for  laterally  homogeneous  earth  models  was  used  to 
fit  both  the  amplitudes  and  travel  times  of  the  observations. 
The  model  obtained  exhibits  a  prominent  low-velocity  zone, 
sharp  discontinuities  of  velocity  at  about  430  and  660  km, 
and  a  discontinuity  in  the  slope  of  the  velocity  at  about 


550  km. 


Long -period  P-Wave 


VI. 3 


Crustal  Structure  beneath  LASA  from 
Spectra  by  W.H.  Bakun 


Abstract 

Long-period  transfer-function  ratios  from  events  in 
South  America,  Fiji-Tonga,  and  Japan  recorded  at  the  LASA 
subarray  centers  are  interpreted  in  terms  of  Haskell- 
Thomson  theory.  The  transfer-function  ratio  data  provide 
a  three-dimensional  model  for  the  crustal  structure  beneath 
LASA.  The  proposed  structure  can  be  characterized  by  two 
trends:  (1)  crustal  thinning  from  the  northeast  to  the 

southwest  across  the  array  and  (2)  a  synclinal  structure 
in  the  southwest  quadrant  of  the  array  with  axis  plunging 
toward  the  northeast. 
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The  response  of  a  crustal  mod?l  of  parallel  plane  layers  of  homo¬ 
geneous  isotropic  material  to  an  incident  ?-;;nve  (K? shell,  was 

utilized  by  Phinney  Wh)  to  define  the  crustal  transfer- function  ratio 
V^)  "  woC^)/U0(-  )•  Wo(*0  and  U0(-0  are  the  spectra  of  the 
vertical  and  horizontal  components  of  displacement  at  the  free  surface 
respectively.  Tp(w)  is  independent  of  the  spectrum  of  the  incident 
P-wave  and  depends  only  on  the  angular  frequency  and  the  apparent 
velocity  c_.  The  position  of  the  peaks  and  troughs  in  Tp(«, )  is 
relatively  insensitive  to  changes  in  phase  velocity  and  sufficiently 
sensitive  to  changes  in  the  crustal  model  parameters  to  be  useful  in 
discriminating  between  competing  crustal  models.  Phinney  proposed  that 
appropriate  models  for  the  crust  and  upper  mantle  beneath  a  recording 
station  could  be  selected  by  matching  )  for  a  suite  of  crustal 
models  to  experimental  transfer- function  ratios  computed  from  the 
Fourier  analysis  of  observed  displacements  or  velocities  of  ?- waves  from 

ISSfSSoSr  Ub0rat0^A°--'-  390  Main  Street, 
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distant  earthquakes.  Ichikawa  and  Bashca  (1965),  Utsu  (1966)5  Sliis  and 
Basham  (1966),  Fernandez  and  Careaga  (1963),  Hasegawa  (1969),  Xurita 

(1969),  and  Bakun  (1970,  1971)  have  utilized  the  transfer-function  ratio 
to  investigate  crustal  structure. 

Churls  (1966)  reported  station  travel- time  residuals  as  large  as 
+1  sec  within  the  200  ton  diameter  of  USA.  Seismic  studies  (Steinhart 
and  Meyer,  1961;  Borcherdt  and  Roller,  19675  Greenfield  and  Sheppard, 
1969;  Glover  and  Alexander,  19 69;  Iyer  at  al..  1970;  Lamer,  1970; 

Mack,  1970;  and  Sngdahl  and  Felix,  1971)  have  confirmed  the  existence 
of  marked  lateral  variations  in  structure  beneath  USA.  While  the 
resultant  crustal  models  differ,  these  studies  agree  that  structure 
beneath  USA  appears  to  be  complex— certainly  not  well  suited  for 

analysis  in  terms  of  the  plane-layer  models  assumed  in  the  Haskell- 
Thomson  theory. 

Bi  this  report,  transfer-function  ratio  data  recorded  at  each  USA 
aubarray  center  are  separated  according  to  the  azimuth  of  P-nave  approach 
in  an  attest  to  isolate  the  effect  of  snail  lateral  areas  beneath  USA. 
As  a  first  approximation,  it  is  assured  that  these  transfer-function 
ratio  data,  indicative  of  sub-US',  structure' beneath  lateral  distance;  of 
approximately  So  to,  can  be  interpreted  in  ten, is  of  parallel  plane¬ 
layered  Haskell-Thonson  models.  The  total  data  imply  a  three-dimensional 
structure  beneath  USA  that  incorporates  important  aspects  of  the 
previously  published  models  for  LASA. 
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DATA 


Transfer- function  ratios  used  in  this  study  were  obtained  from  the 
output  of  the  matched  narrow-band  three- component  long-period  seismometers 
located  at  the  subarray  centers  of  the  A,  C,  D,  2,  and  P  rings  of  LASA. 

The  long-period  seismometers  at  LASA  are  recorded  digitally  with  a 
sampling  intern/ al  of  0.2  seconds.  A  more  complete  description  of  the 
LASA  long-period  system  has  beer,  given  by  Capon  etal.,  (1969).  Ten 
events  (see  Table  1)  with  signal- to-npise  ratios  >2  in  the  0.0?  to  0.1$  Hz 
band  were  selected  from  the  LASA  magnetic  tape  data  library.  (Hoise 
spectra  were  computed  for  each  event  from  a  section  of  record  immediately 
preceding  the  P-wave  onset.)  The  event  population  was  limited  to  events 
in  South  America,  Fiji- Tonga,  and  Japan— the  only  azimuths  from  LASA  for 
which  more  than  one  event  with  an  acceptable  signal- to- noise  ratio-  was 
available.  Typical  data  are  shown  in  Fit  xre  1.  A  30- second  "modified 
Hanning  taper11  window  (Bakun,  1971)  was  applied  to  the  vertical  and  to 
the  resolved  horizontal  longituoinal  time  traces.  The  Fourier  transform 
moduli  of  the  vertical  and  longitudinal  traces  were  then  divided  (e.g.. 
Figure  2)  to  obtain  the  transfer-function  ratio  as  suggested  by  Phinney 
(196U).  Composite  transfer-function  ratios  were  than  formed  for  each 
subarray- azimuth  pair.  The  composite  ratio  is  the  mean  at  each  frequency 
of  the  transfer-function  ratios  available.  The  composite  ratios  for 
each  subarray- azimuth  pair  are  shown  in  Figures  3,  I4,  ?,  and  6.  The 
number  of  events  available  to  define  tne  composite  ratio  for  each 
subarray- azimuth  pair  is  indicated  in  parenthesis  in  Figures  3-6. 


Tho  transfer-function  ratio  samples  structure  near  tho  recording 
site  on  an  azimuth  toward  tho  source  epicenter.  S»  lateral  resolution 
•ia  about  a  crustal  thicknesc  for  a  window  length  of  30  seconds  (Bakun, 
WO)  so  that  crustal  nodal  paraneters  determined  fron  transfer-function 
ratios  at  USA  represent  structure  averaged  over  a  lateral  castar.ce  of 
some  50  Ion  from  the  subarray  center  toward  the  source  epicenter,  ftie 
frequency  of  the  first  (lowest  frequency)  peak  in  the  transfer-function 
ratio  is  diagnostic  of  the  vertical  P-wave  travel  time  through  the 
crust  (Bakun,  1971).  lath  this  in  mind,  the  lowest  freo.uency  peak 
in  the  composite  ratios  shown  in  Figures  3-6  has  been  scored  by  a 
vertical  line.  A  question  mark  adjacent  to  the  vertical  line  indicates 
that  the  composite  ratio  data  for  that  subarroy-azimuth  pair  is  judged 
to  be  of  poorer  quality.  Quality  is  judged  to  be  poor  if:  (a)  only 
one  event  is  available  for  the  subarray-azimuth  pair,  (b)  the  standard 

deviation  of  the  composite  ratio  is  large,  or  (o)  the  resolution  of  the 
lowest  frequency  peak  is  low. 

The  periods  of  tho  lowest  frequency  peak  of  the  composite  ratios 
of  higher  quality  (ratios  not  questioned  in  Figure  e  3-6}  have  been 
Plotted  on  a  map  of  the  USA  array  in  Figure  7.  The  periods  have  been 
Plotted  25  km  (  crustal  thickness)  from  the  subarray  center  on  the 
azimuth  toward  the  source  epicenter.  The  contours  of  the  periods  in 
figure  7  provide  a  three-dimensional  model  for  the  structure  beneath  USA. 
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Piscussira 

Bakun  (1971)  presented  a  relation  between  the  period  of  the  lowest 
frequency  peak  in  Tp(  u))  and  the  vertical  ?-wavc  crustal  travel  tirr.e 
based  on  theoretical  transfer-function  ratios  computed  for  many  crustal 
models.  Using  the  results  of  that  study,  ?-:;avc  times  have  been  assigned 
to  the  contours  in  Figure  7.  Since  crustal  travel  time  increases  rapidly 
with  period  for  periods >  10  seconds,  the  periods  plotted  in  Figure  7 
imply  lateral  differences  of  about  1  second  in  the  vertical  P-wave 
crustal  travel  times  beneath  USA.*  For  example,  assuming  a  constant 
crustal  congressional  velocity  of  6.5  kn/sec,  the  travel  times  in  Figure 
7  imply  an  increase  in  crustal  thickness  from  29  km  in  the  west  to  36  km 
m  the  east,  corresponding  to  an  average  dip  of  the  Koho  of  about  5°, 
plunging  eastward  beneath  LASA. 

The  theoretical  effect  of  a  dipping  crustal  layer  on  the  character 
of  Tp(*j)  is  not  significant  for  dips  less  than  or  equal  to  about  5° 
(Rogers  and  Kisslinger,  1970).  If  velocities  beneath  .LASA  are  greater 
in  the  eastern  sector  than  in  the  west  (Glover  and  Alexander,  1969),  the 
average  dip  on  the  Hoho  beneath  USA  would  tend  to  be  greater  than  5°. 

On  the  other  hand,  all  but  three  of  the  periods  plotted  in  Figure  7  result 
from  updip  approach.  Since  updip  approach  would  increase  the  estimated 
vertical  crustal  travel  times,  or  periods  in  Figure  7,  the  net  effect  of 

the  dip  implied  by  these  data  would  be  to  accentuate  the  lateral  variations 
shown.  • 

The  structure  beneath  USA  shown  in  Figure  7  is  characterized  by 
two  trends:  (1)  crustal  thinning  from  the  northeast  to  the  southwest 
across  the  array  and  (2)  a  synclinal  structure  in  the  southwest  quadrant 
of  the  array  with  axis  plunging  to  the  northeast.  It  is  interesting  to 
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compare  this  structures  obtained  from  transfer- function  ratio  data  with 
previously  published  LASA  models.  The  data  presented  here  tend  ‘bo  con¬ 
firm  Glover  and  Alexander's  (1969)  conclusion  that  the  maximum  lateral 
change  in  structure  beneath  LASA  is  in  a  northeast-southwest  direction, 
with  greater  crustal  thickness  in  the  eastern  sector  of  LASA.  Cn  the 
other  hand,  a  northwest- southeast  profile  through  subarray  AO  would 
result  in  a  two-dimensional  model  similar  in  shape  to  that  proposed  by 
Greenfield  and  Sheppard  (1969  )>  but  displaced  some  SO  km  to  the  south¬ 
east.  Iyer  et  al.  (1970)  and  Larne r  (1970)  proposed  shallower  LISA 
crustal  models  similar  in  shape  to  the  model  of  Greenfield  and  Sheppard. 

The  major  point  of  difference  between  the  model  proposed  here  and 
the  previously  published  crustal  models  for  LASA  lies  in  the  scale  of 
the  structure.  Assuming  reasonable  crustal  velocities,  the  P-wave 
travel  times  of  h.5-5.5  seconds  proposed  here  imply  sub- LASA  crustal 
thicknesses  of  30-1*0  km— significantly  less  than  LASA  crustal  thicknesses 
published  in  previous  studies.  If  the  periods  of  the  lowest  frequency 
peak  of  the  Tp(«J  )  data  were  greater  than  20  seconds,  the  second  (higher 
frequency)  peak  in  Tp(ou )  would  be  erroneously  taken  for  the  desired 
lowest  frequency  pa  ale  in  this  analysis.  Such  a  misidentification  would 
result  in  a  significant  underestimation  of  the  sub- LASA  crustal  times. 


It  should  be  pointed  out  however  that  crustal  thickness  beneath  LASA 
is  a  poorly  determined  parameter.  Seismic  refraction  surveys  (steinhart 
•and  Keyer,  1961;  Borcherdt  and  Roller,  1967)  indicate  crustal  thicknesses 
of  1*7-52  lan  beneath  LASA.  Using  station  travel- time  residuals,  Greenfield 
and  Sheppard  (1969)  proposed  a  model  for  LASA  with  total  crustal  thickness 
of  58  to  70  km.  Larner  (1970)  found  that  structure  similar  in  shape  to 
that  of  Greenfield  and  Sheppard  (1969)  adequately  explained  observed 
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scattering  of  short-period  ?-waves  at  LAS  A..  However,  Larnsr's  analysis 
indicates  an  upper  bound  of  1*8  +p  km  for  the  average  thickness  of  the 
LASA  crust.  Iyer  et  al.  (1970)  found  evidence  for  crustal  thicknesses 
of  1*1  to  55  lot  beneath  LASA. 

It  is  somewhat  surprising  that  a  reasonably  self-consistent  picture 
of  the  complicated  structure  beneath  LASA  can  be  obtained  using  parallel 
plane- layered  theory.  It  should  be  noted  that  only  19  of  the  51  com¬ 
posite  ratios  shown  in  Figures  3-6  were  considered  reliable  enough  for 
use  in  constructing  the  structural  .contours  in  Figure  7.  (identical 
structural  contours—  with  several  inconsistent  data  points— result  if 
all  51  composite  ratios  are  used  in  the  analysis.)  The  19  ratios  used 
here  are  themselves  not  totally  reliable— e.g.,  the  period  of  the  lowest 
frequency  peak  in  Tp(vu)  for  individual  events  in  the  same  subarray- 
azimuth  pair  may  vary.  The  vertical  and  lateral  resolution,  both  on 
the  order  of  a  crustal  thickness,  of  the  transfer- function  ratio  data 
used  here  and  also  the  process  of  forming  composite  ratios  have  the  effect 
of  smoothing  out  more  rapid  variations  in  structure  that  may  be  resolved 
by  short-period  data. 

Scatter  in  the  periods  shown  in  Figure  7  is  most  apparent  between 
the  D  and  2  rings  near  the  El*  subarray,  indicating  rapid  lateral  vari¬ 
ation  in  structure.  It  is  perhaps  only  coincidence  that  this  portion  of 


LASA  lies  atop  Porcupine  Dome,  the  most  significant  geological  feature 
in  the  LASA  area  (King,  1961). 
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Figure  2. 


Figure  3. 
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Figure  5. 
Figure  6. 
Figure  7. 
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VI. 4  Numerical  Simulation  of  Sea-Floor  Spreading  by 
D.J.  Andrews 

Abstract 

Upper  mantle  convection,  including  lithospheric  plate 
motion,  is  simulated  numerically  using  a  two-dimensional 
time -dependent  method  that  allows  large  viscosity  variations. 
The  numerical  operator  developed  for  viscous  flow  is  in  self- 
adjoint  form,  so  that  the  conjugate  gradient  iteration  method 
may  be  used.  Convergence  is  faster  than  with  relaxation  meth¬ 
ods.  The  method  is  applied  to  sea-floor  spreading  with  the 
objective  of  examining  the  driving  mechanism  of  mid-ocean 
ridges.  In  accordance  with  this  objective,  deep  convection 
is  suppressed  in  the  model.  Counterflow  below  the  plates  is 
confined  to  depths  less  than  340  km.  The  model  is  fit  to 
observed  topography  at  four  different  ridge  locations.  It 
is  found  that  for  a  spreading  velocity  of  1.2  cm/yr,  the  ridge 
can  produce  compressive  stress  in  the  lithosphere  out  to  a 
distance  of  1600  km.  For  a  spreading  velocity  of  6  cm/yr 
this  model  is  clearly  excluded,  for  it  requires  excessively 
large  stress  in  the  lithosphere.  Therefore  upwelling  material 
must  cross  the  seismic  discontinuity  at  400  km  depth. 
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A  new  method  is  presented  in  this  paper  for  the  calcu¬ 
lation  of  convection  in  a  fluid  with  strongly  varying 
viscosity.  Time  dependent  heat  flow  and  viscous  mass  flow 
are  simulated  numerically  in  two  space  dimensions.  The 
method  is  designed  for  tectonic  flow  in  the  earth ,  so  that 
inertia  is  neglected.  Therefore  the  Reynolds  number  is 
zero  and  the  Prandtl  number  is  infinite. 

Numerical  calculations  of  steady  state  convection  with 
large  viscosity  variations  have  been  done  by  Torrance  and 
Turcotte  [1971a,  1971b].  They  have  not  included  litho¬ 
spheric  plates  in  their  calculations,  for  extreme  viscosity 
variations  reduce  the  rate  of  convergence  of  the  solution. 
With  the  end  in  view  of  calculating  problems  in  plate 
tectonics,  a  different  approach  to  the  finite  difference 
equations  is  presented  here. 

Equations  to  be  solved 

The  flow  is  assumed  to  be  incompressible,  so  that 
velocity  components  may  be  derived  from  a  stream  function, 

S, 


u  = 


ds 

3y 


(1) 
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Components  of  the  stress  deviator  tensor  are  assumed  to  be 
determined  by  Newtonian  viscosity  in  two  dimensional  flow. 
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The  viscosity,  n,  need  not  be  constant.  Pressure  (mean 
stress)  may  assume  any  value  needed  to  produce  stress 
equilibrium.  In  the  absence  of  inertia  each  component  of 
force  density  vanishes 
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We  adopt  the  Boussinesq  approximation,  that  only  the 
variation  of  density  with  temperature  need  be  considered, 
and  only  in  the  body  force  term.  A  constant  thermal 
expansivity  is  assumed 

P  ■  PQ  (1  -  «T) 
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(7) 
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although  this  relation  could  be  generalized  with  no 

difficulty.  Overburden  pressure  due  to  a  constant  density 
is  subtracted  from  pressure , 

p  =  p  -  pQgy 

Then  equations  (5)  and  (6)  become 


where  all  terms  have  about  the  same  order  of  magnitude. 
Pressure  may  be  eliminated  by  taking  the  curl  of  force 
density.  The  single  equation  that  results,  expressed  in 
terms  of  the  stream  function,  is 


(ID 


For  a  fixed  temperature  field,  mass  flow  will  be 
steady.  Time  dependence  of  the  solution  arises  from  time* 
dependence  of  temperature.  At  each  instant  of  time  the 
finite  difference  approximation  of  this  elliptic  equation 
is  solved  by  an  iterative  procedure. 

Relaxation  methods  for  solution  of  systems  of  linear 
equations  can  be  put  upon  a  rigorous  basis  only  for 
symmetric  matrices  [Varga ,  1962] ,  although  such  methods  have 
been  used  more  generally.  In  the  hope  of  handling  extreme 
viscosity  variations,  the  present  worV.  will  be  kept  as 
well-founded  as  possible. 

The  differential  operator  appearing  in  equation  (11) 
is  self-adjoint  and  can  be  derived  from  a  variational 
principle  of  minimum  viscous  dissipation.  If  a  finite 
difference  analogue  of  this  operator  is  formulated  from  a 
variational  principle,  then  one  has  increased  confidence 
that  an  iterative  relaxation  procedure  will  converge.  Also 
it  will  be  seen  that  the  finite  difference  operator  so 
derived  is  a  self-adjoint  (symmetric)  matrix. 

The  single  fourth  order  equation  for  stream  function 
might  be  replaced  by  a  set  of  coupled  equations  of  lower 
order  in  more  than  one  dependent  variable.  For  instance, 
velocity  components  might  be  used  rather  than  stream  function 
as  dependent  variables.  One  might  expect  a  variational 
principle  to  hold  with  incompressibility  imposed  as  a 
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subsidiary  condition.  However,  a  divergence  expression 
cannot  be  used  as  a  subsidiary  condition  in  a  variational 
problem  [Courant  and  Hilbert,  1953].  Alternatively,  a 
coupled  pair  of  equations  could  be  used  for  vorticity 
and  stream  function.  One  would  like  to  use  a  relaxation 
procedure  for  each  equation  of  the  pair  [Burggraf ,  1966], 
but  in  the  case  of  variable  viscosity  the  differential 
operators  of  the  separate  equations  cannot  be  derived 
from  variational  principles.  Vorticity  has  lost  its  use¬ 
fulness  as  an  auxilliary  variable.  Therefore  it  is 
desirable  to  work  with  the  single  equation  of  fourth  order 
for  the  stream  function. 

However,  there  is  a  disadvantage  to  using  a  relaxation 
procedure  with  a  fourth  order  equation.  By  dimensional 
arguments  one  would  expect  the  number  of  iterations 
required  for  each  component  of  the  error  to  decay  to  be 
proportional  to  the  fourth  power  of  its  wavelength.  Then 
che  number  of  iterations  required  would  go  as  the  fourth 
power  of  the  number  of  finite  difference  zones  in  one 
space  dimension,  rather  than  the  square  as  for  second  order 
equations.  Fortunately,  with  the  symmetric  matrix  equation 
formulated  here,  one  can  use  the  conjugate  gradient  method, 
that  ideally  gives  an  exact  solution  after  a  number  of 
iterations  proportional  to  the  square  of  the  number  of 
zones  in  one  dimension. 
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The  aquation  used  for  evolution  of  the  temperature 
field  in  time  is 
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(12) 


where  K  is  thermal  conductivity  and  H  is  the  rate  of 
heat  generation  per  unit  volume .  The  terms  on  the  right 
side  account  for  advection  of  heat  by  mass  transport, 
thermal  conduction,  viscous  dissipation,  adiabatic 
temperature  change  with  depth,  and  radioactive  heat 
generation . 

Finite  Difference  Equations 

In  describing  the  finite  difference  equations,  inter¬ 
sections  of  grid  lines  will  be  called  grid  points  and  mesh 
spaces  between  grid  lines  will  be  called  zones  or  zone 
interiors.  Grid  lines  are  parallel  to  the  x  and  y  coordinate 
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axes.  The  stream  function  is  defined  at  grid  points.  Then 
the  stress  deviator 

8xx  "  2n 


is  most  naturally  defined  by  a  finite  difference  expression 
centered  in  a  zone  interior 


(sxx'i-i,j-i  "  2"i-i. j-i  (Si.3  ' 
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where  subscripts  are  column  number  and  row  number  in  the 
grid. 

The  stress  component 
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is  most  naturally  defined  by  a  finite  difference  expression 
centered  at  a  grid  point.  For  the  case  of  square  zones  of 
uniform  size  the  expression  is 

<s*y>i,j  -  ‘  +  Si-^1 
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Of  course ,  strain  rates  are  given  by  the  same  expressions 
divided  by  2r\ . 

The  finite  difference  expressions  for  strain  rates  may 
be  multiplied  by  finite  difference  expressions  for  stress 
deviators  and  summed  to  get  an  expression  for  shear 
dissipation.  In  the  continuous  case  the  differential  operator 
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may  be  found  by  minimizing  the  volume  integral  of  shear 
dissipation.  We  will  form  a  finite  difference  analogue 
of  this  differential  operator  by  minimizing  our  finite 
difference  expression  for  shear  dissipation.  The  value 
of  the  stream  function  at  a  particular  grid  point  enters 
into  expressions  for  shear  dissipation  at  the  four 
neighboring  grid  points  and  four  neighboring  zone  interiors. 
The  sum  of  these  terms  is  minimized  with  respect  to  the 
stream  function  at  the  central  grid  point  when  the 
following  equation  is  satisfied. 


4na  (SNE  -  SN  +  s0  -  SE) 

‘  4VSN  "  SNW  +  SM  '  S0> 

+  4VS0  -  +  ssw  -  ss) 
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nE(SNE  +  SSE  "  SEE  ‘  s0> 
nK(SUW  +  SSW  ~  so  ~  sww)  ”  0 


where  subscripts  indicate  pcsitions  shown  in  Figure  1.  it 
is  assumed  that  the  zones  are  square  and  uniform  in  size. 

More  generally  there  will  be  a  term  on  the  right-hand  side 
arising  from  buoyant  body  force. 

This  finite  difference  analogue  of  the  differential 
operator  not  only  satisfies  a  finite  difference  variational 
Principle  precisely,  but  also  gives  a  solution  that  precisely 
satisfies  analogues  of  the  stress  equilibrium  equations  (9) 
and  (10)  with  Stress  components  defined  by  equations  (13)  and 
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(14).  Viscosity  may  be  defined  arbitrarily  at  each  grid 
point  and  zone  interior.  Moreover,  the  self-adjoint 
property  of  the  differential  operator  is  preserved,  for  the 
matrix  relating  values  of  stream  function  at  different 
grid  points  is  symmetric.  This  may  be  seen  more  readily 
if  the  expression  is  rearranged  as  follows, 

[4(na  +  nb  +  nc  +  na>  +  nE  +  nN  +  nw  +  nsis0 

-  «Oa  +  Vse  -  4(na  +  nb)SN 

-  4(nb  +  ti0)sw  -  4 (nc  +  na)ss 

+  (4na  "  nE  "  nN,SNE  +  <4nb  *  r,N  "  nM)SNW 

+  (4n0  -  1W  -  Is)ssw  "  (4na  ■  ns  '  VSSE 

+  nESEE  +  ^N^NN  +  nWSWW  +  I1SSSS 

-  -  <eo9C»0  (Ta  -  Tb  -  Tc  +  Ta)  <4x)*/2 

Here  the  buoyant  term  has  been  written  on  the  right  hand 
side.  A  system  of  such  equations  is  to  be  solved,  one 
equation  for  each  grid  point,  with  the  diagram  of  Figure  1 
shifted  so  that  it  is  centered  on  each  point.  The  matrix 
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element  linking  S^,  to  Sq  is  -4  { n ^  +  Tia) .  If  the  diagram 

is  shifted  to  the  right  one  -.one,  one  can  see  that  the  matrix 

element  linking  the  same  two  values  (now  called  SQ  and  Sw) 

is  -4  (ru  +  ti  )  ,  which  is  the  same.  One  may  examine  the 
d  c 

other  terms  in  the  same  way  to  show  that  the  matrix  is 
symmetric . 

We  have  chosen  to  define  temperature  in  zone  interiors, 
that  is,  at  points  staggered  with  respect  to  points  where 
the  stream  function  is  defined.  In  formulating  a  finite 
difference  analogue  of  equation  (12)  the  advection  terms 
must  be  considered  carefully.  In  a  mantle  convection 
problem  the  advection  terms  will  be  several  times  larger 
than  the  conduction  terms,  so  it  is  important  that  no 
anomalous  diffusion  is  introduced  by  the  finite  difference 
approximation . 

If  a  system  of  implicit  finite  difference  equations 
were  used  for  temperature,  the  advection  terms  would 
contribute  antisymmetric  components  to  the  matrix  operator. 
Alternating  direction  implicit  (ADI)  methods,  as  well  as 
relaxation  methods,  are  well-founded  only  for  symmetric 
matrices.  Although  numerical  methods  are  commonly  used 
successfully  in  situations  that  cannot  be  rigorously 
defended,  a  more  conservative  approach  will  be  used  here. 
Explicit  finite  difference  expressions  will  be  used  for  the 
advection  terms.  To  insure  stability  the  time  step  must  be 


less  than  the  time  for  a  fluid  particle  to  cross  one  zone. 

For  a  zone  size  of  10  or  20  kilometers  and  a  fluid  velocity 
of  a  few  centimeters  per  year,  the  stability  criterion 
imposed  by  explicit  conduction  terms  is  no  more  stringent 
than  this  advection  criterion.  Therefore  the  entire  finite 
difference  approximation  to  equation  (12)  will  be  explicit. 

In  the  following  equations  the  superscript  o  for 
old  indicates  a  value  from  the  previous  time  step,  the  super¬ 
script  n  for  new  indicates  a  value  updated  to  the  current 
step  in  time,  and  h  indicates  a  value  half-way  through 
the  time  step.  To  calculate  the  advection  terms  temperatures 
are  first  defined  on  faces  between  zones  half-way  through 
the  time  step, 

TAh  =  (T0°  +  T^J/2  +  uA°  (Tj®  -  T0°)At/(2Ax)  (17) 

Tfih  =  (Tq°  +  T2°)/2  +  vb°(T2°  -  Tg°) At/ (2Ay)  (18) 

where  subscripts  indicate  positions  in  Figure  2.  Then  an 
updated  temperature  due  to  advection  is  found  by  an  equation 
in  conservation  form 

T0n  -  T0°  -  (uA°TAh  -  uc°Tch)At/i* 

-  (uB°TBh  ‘  UD°TDh)4t/AX 


(19) 
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This  procedure  is  accurate  to  second  order.  For  a  uniform 
velocity  field  there  is  second  order  dispersion  and  fourth 
order  diffusion  of  temperature,  as  for  the  Lax-Wendroff 
procedure . 

Finite  difference  approximations  to  the  remaining 
terms  of  equation  (12)  are  not  novel  and  will  not  be  listed 
here. 

After  temperature  is  advanced  a  step  in  time,  new 
values  of  viscosity  and  updated  buoyant  terms  can  be 
evaluated.  Then  an  iteration  must  be  performed  to  find 
the  stream  function  appropriate  to  the  new  temperature 
field.  A  few  relaxation  iterations  suffice  to  correct 
short  wavelength  components .  After  a  long  wavelength  error 
has  accumulated,  the  complete  conjugate  gradient  iterative 
procedure  should  be  performed. 

The  Method  of  Conjugate  Gradients 

The  solution  of  the  system  of  equations  for  the  stream 
function  over  a  two-dimensional  grid  will  be  discussed  in 
this  section  as  if  the  unknown  values  were  re-ordered  into 
a  vector  s  of  length  N.  Then  the  system  of  equations  to 
be  solved  is  represented  by  the  matrix  equation 

K  s  =  b 

(PRINTER:  Single  underline  in  bold  face;  Double  underline  in 
bold  face  sans  serif) 
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The  vector  on  the  right-hand  side  consists  of  the  buoyant 
terms  and  terms  containing  fixed  boundary  values.  Boundary 
values  proportional  to  unknown  values  inside  the  grid  are 
treated  as  modifications  to  some  components  of  the  matrix 
K.  The  matrix  K  for  the  viscous  flow  equations  formulated 
in  the  previous  section  is  symmetric ,  and  the  boundary 
conditions  that  have  been  used  preserve  the  symmetry.  Also, 
it  is  physically  reasonable  that  K  is  positive  definite, 
although  no  proof  has  been  found.  It  is  not  necessary  to 
store  the  entire  matrix  K,  since  all  components  relating  to 
a  particular  component  of  s  can  be  formed  easily  from 

viscosity  values  in  the  neighborhood. 

Bennett  and  Cooley  [1969]  devised  a  method  of  solving 
such  a  system  of  equations  that  does  not  require  that  an 
NXN  matrix  be  formed  and  stored.  Their  method  is  based 
on  earlier  ideas  of  Hestenes  and  Stiefel  [1952]  and 
Lanczos  [1952] .  It  is  assumed  that  K  is  symmetric  and 
positive  definite.  Let  s1  be  the  approximation  to  the 
solution  in  the  ith  iteration  and  let  r1  be  its 
residual, 


Also,  an  auxilliary  vector^  z1  appears  in  the  algorithm. 
The  initial  step  of  the  algorithm  is 


Then  the  following  sequence  of  equations  is  executed 
iteratively, 


g1  »  -  (riT  z1)/^11  K  z1) 


i+1  ,  i  i  i 
s  *  s  -  g  z 


„i+l  i  .  1  v  _i 

r  £  +  9  —  — 


f1  -  (ri+1T  ri+1)/(riT  rL) 


zi+1  -  ri+1  +  f1  z1 


In  each  iteration  four  vectors  must  be  stored,  s  ,  r  , 

•  • 
z1,  and  £  z1. 

The  reader  is  referred  to  Bennett  and  Cooley  [1969] 
for  a  proof  of  the  algorithm.  Residuals  from  different 
iterations  are  orthogonal 


riT  rj  .  o  ,  i  +  j 


and  the  vectors 


are  conjugate  with  respect  to  K  , 
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More  importantly,  the  vectors  z^  are  linearly  independent, 
so  that  the  successive  corrections  to  the  stream  function, 
gizi,  span  the  entire  N  dimensional  space,  after  N 
iterations.  Also  the  coefficient  of  each  correction,  g1 , 
given  by  equation  (23)  ,  minimizes  ei+lT  £  e1+1 ,  where 
e1+1  is  the  error  s  -  si+1.  Therefore,  if  there  were 
no  round-off  error,  an  exact  solution  would  be  reached  in 
N  iterations  or  less. 

To  reduce  the  affect  of  round-off  error  it  is 
desirable  that  the  error  of  the  initial  approximation  have 
negligible  components  with  large  eigenvalues  (short 
wavelengths) .  This  can  be  accomplished  with  a  few  relax¬ 
ation  iterations  before  beginning  the  conjugate  gradient 
procedure . 

In  our  viscous  flow  calculations  it  has  been  found  that 
the  calculated  product  ziT  g  z1  is  always  positive, 
giving  support  to  our  conjecture  that  our  g  is  positive 
definite.  Long  wavelength  components  of  the  solution 
converge  faster  with  the  conjugate  gradient  method  than  with 
relaxation  methods.  On  an  IBM  360  with  six  significant 
figures  convergence  has  been  satisfactory  with  viscosity 
varying  by  a  factor  of  1000.  With  double  precision 
arithmetic  larger  viscosity  contrasts  have  been  used  with 
no  difficulty. 
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Mid-Ocean  Ridge  Calculations 

The  general  calculational  method  of  the  first  part  of 
this  paper  is  used  here  to  investigate  a  single  aspect  of 
plate  tectonics,  namely  the  effectiveness  of  a  mid-ocean 
ridge  as  a  driving  mechanism.  General  considerations 
of  plate  motion  driven  by  thermally  induced  density 
instabilities  have  been  given  by  McKenzie  [1969] .  In 
discussions  of  driving  mechanisms,  attention  is  commonly 
focused  either  on  the  plates  themselves,  or  on  convection 
calls  below  the  plates. 

In  the  first  class  of  models,  lithospheric  plates 
themselves  are  recognized  as  an  important  part  of  the  total 
mass  circulation  [Elsasser.  1971] ,  [Jacoby.  1970] .  The 
density  contrast  between  lithosphere  and  asthenosphere  is 
a  driving  mechanism.  Plates  are  pushed  by  ridges  and  pulled 
by  down-going  slabs.  Plates  are  rigid  enough  compared  to 
the  asthenosphere  to  act  as  stress  guides.  They  may  be 
decoupled  by  the  asthenosphere  from  the  rest  of  the  mantle. 

In  the  second  class  of  models  the  mantle  below  the 
lithosphere  is  considered  to  be  convecting,  and  plates  ride 
passively  on  the  tops  of  convection  cells  [Torrance  and 
Turcotte,  1971a,  1971b]. 

These  two  classes  should  not  be  considered  antithetical, 
for  they  both  may  be  important  aspects  of  the  complete 
picture.  The  model  considered  here  belongs  to  the  first 


class  with  attention  focused  on  the  lithosphere.  However, 
it  must  be  recognized  that  convection  is  taking  place 

through  much  of  the  mantle. 

A  Newtonian  fluid  of  constant  viscosity  is  unstable 
against  convection  if  the  Rayleigh  number  exceeds  a 
critical  value.  Recall  that  the  Rayleigh  number  increases 
with  larger  heat  flow,  smaller  conductivity,  smaller 
viscosity,  and  greater  depth  of  convection  cells.  Recent 
determinations  of  viscosity  by  Cathles  11971]  and  of 
conductivity  by  Schatz  [1971]  strongly  support  mantle-wide 
convection.  Cathles  concludes  from  a  comprehensive 
analysis  of  glacial  uplift  data  that  viscosity  does  not 
exceed  1022  poise  anywhere  in  the  mantle.  With  Cathles' 
viscosity,  and  any  reasonable  values  for  other  physical 
parameters,  the  critical  Rayleigh  number  is  exceeded  for 
the  mantle  as  a  whole. 

Schatz  [1971]  has  determined  that  the  radiative 
portion  of  thermal  conductivity  of  olivine  is  not  as 
important  at  high  temperature  as  had  previously  been 
thought.  He  concludes  that  conductivity  in  the  mantle  is 
roughly  constant  as  a  function  of  depth.  If  heat  were 
transported  by  conduction  alone  in  the  mantle,  Schatz 's 
value  of  conductivity  would  require  temperatures  far  above 
the  solidus  for  most  of  the  mantle.  Furthermore,  with 
Cathles'  viscosity  and  Schatz 's  conductivity,  convection 
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may  take  place  with  cells  only  140  km  from  top  to  bottom. 
Therefore  convection  is  certainly  taking  place. 

In  this  work  we  will  consider  plate  movement  away 
from  the  ridge  axis  and  a  counterflow  below  the  plate 
confined  to  depths  less  than  340  km.  The  objective  is  to 
consider  this  relatively  shallow  flow  and  counterflow 
independently  of  the  deeper  circulation  that  is  certainly 
taking  place,  and  to  establish  the  conditions  for  which 
this  model  is  reasonable. 

In  order  to  separate  effects,  it  is  necessary  to 
artificially  inhibit  convection  below  the  lithosphere  in 
our  model.  The  initial  geotherm  used  should  be  stable 
against  convection  in  the  depth,  considered.  To  achieve 
this  objective  an  artificially  high  conductivity  of  0.1 
cal/cm-sec-deg  is  used  for  temperatures  above  1275  °C. 

In  this  way  the  heat  flux  from  the  earth's  interior  is 
obtained  by  conduction  alone  with  a  small  thermal  gradient 
(1  deg/km) ,  about  the  magnitude  to  be  expected  in  the 
case  of  convection.  A  realistic  conductivity  is  used  for 
the  lithosphere,  0.0075  cal/cm-sec-deg. 

—6  2 

A  normal  surface  heat  flux  of  1.1  x  10  cal/cm  -sec 

-14 

is  assumed  and  radioactive  heat  generation  of  1  x  10 
cal/cm  -sec  is  adopted  [Sclater  and  Francheteau,  1970] . 

The  normal  geotherm  corresponding  to  these  parameters  is 
shown  in  Figure  3a.  It  is  assumed  that  radioactivity  is 


uniform,  for  no  chemical  differentiation  is  accounted  for  in 
this  model. 

The  viscosity  model  used  is  shown  in  Figure  lb  for 
the  normal  geotherm.  Values  of  viscosity  for  the  astheno- 
sphere  and  mesosphere  are  taken  from  Cathles  [1971]  . 

The  transition  from  lithosphere  to  asthenosphere,  normally 
at  80  km,  [Press  and  Kanamori,  1970]  is  prescribed  to  take 
place  at  1130  #C.  The  transition  from  asthenosphere  to 
mesosphere  is  fixed  at  160  km  depth.  The  lithosphere  is 
treated  as  a  viscous  fluid  with  viscosity  much  larger  than 
the  rest  of  the  mantle.  The  actual  value  of  viscosity  used 
in  the  lithosphere  is  not  important  in  this  calculation, 
if  it  is  large  enough  to  ensure  nearly  rigid  motion. 

A  closed  circulation  path  is  not  considered  in  this 
work.  An  artificial  lateral  boundary  is  used,  on  which 
horizontal  input  and  output  velocities  are  prescribed. 

From  the  calculation,  stress  at  the  boundary  is  found. 
Although  it  is  artificial  to  separate  this  region  from  the 
rest  of  the  world.,  stress  at  the  boundary  gives  a  quanti¬ 
tative  indication  of  their  mutual  interaction  for  the 
assumed  spreading  velocity. 

The  prescribed  velocity  profile  at  the  lateral 
boundary  is  shown  in  Figure  3c.  This  profile  was  chosen  to 
minimize  dissipation  for  one-dimensional  shear  flow.  This 
profile  is  scaled  to  get  whatever  plate  velocity  io  desired 
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From  this  velocity  profile  shear  drag  at  the  base  of  the 

lithosphere  can  be  found.  It  is  proportional  to  velocity 

and  is  9  bars  for  1.2  cm/yr  half  spreading  rate.  The  two- 

dimensional  calculations  show  nearly  horizontal  flow  and 

smooth  variation  of  stress  near  the  lateral  boundary, 

thus  supporting  the  treatment  of  flow  through  this  boundary. 

A  constant  thermal  expansivity,  a  ■  4  x  10  5  deg  1, 

is  used  [Sclater  and  Francheteau,  1970;  Sleep,  1969]. 

This  value  affects  flow  through  buoyancy  and  affects 

topography  by  thermal  expansion  of  the  lithosphere.  Other 

physical  parameters  are  heat  capacity,  Cp  =  0.311  cal/gm-dfeg, 

density,  pQ  -  3.4  gm/cm3,  and  acceleration  of  gravity, 

2 

g  ■  990  cm/sec  . 

Each  of  the  calculations  was  done  with  finite  difference 
zone  size  of  20  km,  and  with  the  calculational  region 
extending  to  a  depth  of  340  km  and  a  distance  of  1000  km 
from  the  ridge  axis.  The  ridge  axis  was  treated  as  a  plane 
of  reflection. 

Results 

Streamlines  and  isotherms  are  shown  in  Figure  4  for 
the  case  of  1.2  cm/yr  plate  velocity  run  to  steady  state. 
Streamlines  within  the  lithosphere  are  quite  close  to  being 
straight  and  parallel.  The  velocity  of  upwelling  at  the 
ridge  axis  is  1.6  cm/yr,  greater  than  the  plate  velocity 


112 


in  this  case.  The  flow  pattern  is  not  similar  for  different 
spreading  rates.  For  a  plate  velocity  of  6.0  cm/yr  the 
pwelling  velocity  is  2.6  cm/yr.  Shear  heating  far  from 
the  ridge  in  the  case  of  1.2  cm/yr  plate  velocity  con¬ 
tributes  only  1.5%  to  the  normal  surface  heat  flux.  Shear 
heating,  of  course,  is  proportional  to  the  square  of  plate 
velocity  for  moderate  velocities,  and  then  increases  less 
rapidly  since  the  increased  heat  flow  thins  the  lithosphere. 

The  time  dependent  calculations  done  represent 
spreading  histories  at  four  different  ridge  locations.  A 
calculation  to  steady  state  was  run  for  1.2  cm/yr  to 
represent  the  Mid-Atlantic  ridge  at  32  N.  This  spreading 
rate  was  averaged  from  Pitman  at  al.  [1971] .  A  constant 
velocity  of  1.8  ciVyr  was  used  to  represent  the  Mid-Atlantic 
ridge  at  28  S,  [Dickson,  1968;  Maxwell ,  1969].  A  variable 
spreading  rate  was  used  for  the  Reyk janes  ridge  of  1.2  cm/yr 
for  0  to  9  m.y.b.p. ,  0.4  cm/yr  for  9  to  37  m.y.b.p.,  1.8 
cm/yr  for  37  to  50  m.y.b.p.,  and  1.2  cm/yr  earlier  [Pitman 
al* 1  1971] .  The  East  Pacific  Rise  at  43  S  is  represented 
by  spreading  history  of  6  cn\/yr  for  0  to  9  m.y.b.p.  and 
1.8  cin/yr  before  that  [Pitman,  1968].  According  to  Pitman's 
magnetic  anomaly  map  the  ridge  is  spreading  asymmetrically. 
The  spreading  history  used  here  applies  to  the  east  flank. 

Surface  heat  flow  found  in  each  of  these  calculations 
closely  agrees  at  corresponding  ages  of  ocean  floor.  Since 
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observed  heat  flow  data  have  large  scatter,  it  is  reasonable 
to  make  one  plot  applicable  to  all  locations  by  using  age 
as  the  independent  variable.  In  Figure  5  calculated  heat 
flow  is  compared  with  data  from  the  Pacific  tsclater  and 
Francheteau ,  1970).  Calculated  heat  flow  values  for  the 
two  steady  state  cases  ere  identical  at  corresponding  ages 
„d  are  shown  as  Curve  A.  Calculated  heat  flow  for  the 
Reykjanes  and  East  Pacific  ridges,  Curves  B  and  C,  deviate 
from  the  steady  state  case  much  less  than  the  deviation 
of  observed  data.  The  error  bar.  indicate  standard  devi¬ 
ations  of  the  observations.  Since  measurements  are  made 
preferentially  in  sediment  ponds,  the  data  may  be  as  much 
a.  50%  too  low  tsclater  and  Francheteau,  1970).  If  there 
is  water  circulation  through  fractures  in  the  ocean  floor, 
a  significant  part  of  the  total  heat  flow  may  not  be 
measured,  especially  in  rough  topography  near  the  ridge 
crest  r Lister,  1971).  Therefore  it  is  not  unreasonable 
that  the  theoretical  prediction  is  at  the  upper  side  of  the 

scatter  of  data. 

in  the  viscous  flow  calculation  the  upper  boundary  is 
a  rigid  horizontal  plane  on  which  free  sliding  is  allowed. 
Prom  the  calculated  stream  function,  pressure  and  viscous 
stress  deviators  are  known  throughout  the  region.  The 
vertical  stress  component  (pressure  plus  deviator)  is  not 
sero  at  the  upper  surface.  From  this  component  one  obtains 
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mass  per  unit  area  lying  above  this  reference  plane,  and 
hence  elevecion  of  the  ocean  floor. 

In  Figures  6  through  9  calculated  results  for  the 
four  ridge  locations  are  compared  to  observed  data.  In 
each  case  an  observed  topography  is  shown  [Talwani ,  1970]  , 
together  with  calculated  topography  resulting  from  viscous 
stress.  Also  shown  is  topography  calculated  from  temper¬ 
ature  assuming  isostatic  compensation  at  100  km  depth. 

Comparison  with  free  air  gravity  measurements  is  not 
shown,  since  the  prediction  is  sensitive  to  temperature 
and  stress  at  the  artificial  lower  boundary.  The  fact 
that  observed  mean  free  air  anomalies  are  near  zero  may  be 
a  sensitive  test  for  later  work  on  the  entire  mantle 

•  4 

circulation.  Predicted  Bouguer  anomalies,  dependent 
primarily  on  topography,  are  compared  to  observed  for 
the  Reyk janes  [Talwani  et  al.  1971]  and  North  Atlantic 
ridges  [Talwani  et  al. ,  1965] . 

In  the  1.2  and  1.8  cii^/yr  cases  topography  calculated 
from  viscous  stress  and  isostatic  topography  differ  only 
slightly.  In  this  case  the  cooling  slab  model  of  McKenzie 
[1967]  and  Sleep  [1969]  is  valid  for  topography  and  heat 
flow.  As  in  that  model,  one  can  adjust  conductivity  and 
expansivity  to  fit  topography,  but  it  is  not  possible  to 
fit  both  topography  and  mean  heat  flow  with  one  set  of 
parameters  [Sleep,  1969;  Sclater  and  Francheteau.  1970]. 
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In  the  case  of  the  East  Pacific  Rise  the  calculated 
topography  is  far  from  being  isostatic. 

The  effectiveness  of  a  mid-ocean  ridge  as  a  driving 

mechanism  is  demonstrated  in  Figure  10 ,  which  shows  the 

force  transmitted  by  the  lithosphere,  acting  as  a  stress 

guide.  The  horizontal  stress  deviator  component  is 

integrated  through  the  thickness  of  the  lithosphere,  and 

the  result  is  expressed  in  terms  of  the  equivalent  stress 

for  a  thickness  of  60  km.  Beyond  1000  km  reasonable 

extrapolations  are  shown.  One  might  ask  why  the  viscous 

contribution  to  mean  stress,  the  nonisostatic  pressure, 

is  not  taken  into  account  in  the  Figure.  In  all  cases 

calculated  the  work  done  by  nonisostatic  pressure  at  the 

lateral  boundary  was  nearly  the  same  as  the  work  done  by 

the  stress  deviator  there.  Therefore  the  stress  deviator 

is  a. good  indicator  of  driving  mechanism. 

% 

Drag  on  the  base  of  the  lithosphere  requires  that 
horizontal  stress  become  more  tensile  with  increasing 
distance  from  the  ridge  axis.  This  explains  the  positive 
slope  of  force  with  respect  to  distance  far  from  the  ridge. 
Closer  in,  where  the  lithosphere  is  gradually  thickening 
compressive  stress  is  built  up  due  to  the  buoyancy  of 
underlying  material.  For  the  slow  moving  ridges  this 
mechanism  is  an  important  part  of  the  overall  driving 
mechanism.  At  a  velocity  of  1.2  cm/yr,  the  ridge  can 
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produce  compressive  stress  in  the  lithosphere  out  to  a 
distance  of  1600  km. 

For  the  East  Pacific  Rise,  however,  the  horizontal 
stress  deviator  in  the  lithosphere  is  tensile  at  all 
distances,  and  reaches  unrealistically  high  values  within 
1000  km.  This  tensile  stress  would  be  high  enough  to 
produce  continual  earthquakes  in  the  interior  of  the  pla^-a. 
Such  earthquakes,  of  course,  are  rare.  Mendiguren  [1971] 
analyzed  the  source  mechanism  of  a  shock  from  the  middle 
of  the  Nazca  plate,  and  found  horizontal  compressive  stress, 
in  contradiction  to  the  prediction  of  this  model  for  fast 
spreading . 

Therefore  the  counterflow  for  fast  moving  plates 
cannot  be  confined  to  relatively  shallow  depths,  as 
postulated  here.  The  driving  mechanism  in  these  cases 
certainly  must  involve  upwelling  from  greater  depths,  and 
possibly  a  convection  cell  of  considerable  depth  in  the 
mantle.  The  asthenosphere  is  not  weak  enough  to  decouple 
plates  from  flow  lower  in  the  mantle. 

Acknowledgements 

The  author  is  grateful  to  John  Minear  for  calling 
his  attention  to  the  method  of  conjugate  gradients,  to 
Norman  Sleep  and  Larry  Cathles  for  many  stimulating  dis¬ 
cussions  and  to  Jean-Claude  de  Bremaecker  for  our  mutual 


117 


exchange  of  ideas  in  this  area.  Judy  Andrews  helped 
with  the  library  search  and  preparation  of  figures. 

This  research  was  supported  by  a  grant  to  Professor 
M.  Nafi  Toksdz  by  the  Advanced  Research  Projects  Agency 
and  monitored  by  the  Air  Force  Office  of  Scientific 
Research  under  contract  AF49  (638) -1763 . 


REFERENCES 


Bennett,  J.M.  and  P.C.  Cooley,  Linear  analysis  of 

structures  by  the  method  of  conjugate  gradients , 

41  pp. ,  Basser  Computing  Department,  University  of 
Sydney,  Technical  Report  51.  1969. 

Burggraf,  O.R.,  Analytical  and  numerical  studies  of  steady 
separated  flows,  J.  Flu.  Mech..  2A,  113,  1966. 

Cathles,  Lawrence,  The  viscosity  of  the  earth's  mantle, 
PH.D.  thesis,  Princeton  University,  1971. 

Courant,  R.  and  D.  Hilbert,  Methods  of  Mathematical 

Physics,  Vol.  I,  561  pp. ,  Interscience  ,  New  York, 
1953. 

Dickson,  G. ,  E.  Pitman,  and  J.  Heirtzler,  Magnetic  anomalie 
in  the  south  Atlantic  and  ocean-floor  spreading, 

J .  Geophys .  Res . ,  73,  2087,  1968. 

Elsasser ,  W.M.,  Sea  floor  spreading  as  thermal  convection, 
J.  Geophys.  Res..  76*  1101,  1971. 

Hestenes,  M.,  and  E.  Stiefel,  Methods  of  conjugate  gradient 
for  solving  linear  systems,  Journal  Res.  NBS.  49, 

409,  1952. 

Jacoby,  W.R. ,  Instability  in  the  upper  mantle  and  global 
plate  movements,  J .  Geophys .  Res . ,  75,  5655,  1970. 

Lanczos,  C.,  Solution  of  systems  of  linear  equations  by 

minimized  iterations,  Journal  Res.  NBS.  49,  33,  1952. 


Lister,  C.R.B.,  On  the  thermal  balance  of  a  mid-ocean  ridge, 
(abstract),  Trans.  Am.  Geophys.  Union.  52 ,  923,  1971. 

McKenzie,  D.P.,  Some  remarks  on  heat  flow  and  gravity 
anomalies,  J .  Geophys .  Res . .  72,  6261,  1967. 

McKenzie,  D.P.,  Speculations  on  the  consequences  and 

causes  of  plate  motions,  Geophys.  J.  R.  astr.  Soc. . 

18,  1,  1969. 

Maxwell,  E.E.,  R.P .  Von  Herzen,  K.J.  Hsu,  J.E.  Andrews, 

T.  Saito,  S.F.  Percival,  Jr.,  E.D.  Metow,  and  R.E. 
Boyce,  Recent  deep  sea  drilling  results  from  the 
South  Atlantic,  Science,  168,  1047,  1970. 

Mendiguren,  J.A. ,  Focal  mechanism  of  a  shock  in  the  middle 
of  the  Nazca  plate,  J.  Geophys.  Res..  76,  3861,  1971. 

Pitman,  w.c.  Ill,  E.M.  Herron,  and  J.R.  Heirtzler,  Magnetic 
anomalies  in  the  South  Pacific  Ocean  and  sea-floor 
spreading,  J .  Geophys .  Res . ,  73 ,  2069,  1968. 

Pitman,  W.C.  Ill,  M.  Talwani,  and  J.R.  Heirtzler,  Age  of 

the  North  Atlantic  Ocean  from  magnetic  anomalies.  Earth 
and  Planet.  Sci.  Letters.  11,  195,  1971. 

Press,  F.,  and  H.  Kanamori,  How  thick  is  the  lithosphere, 
(abstract),  Trans.  Am.  Geophys.  Union.  51,  363,  1970. 

Schatz,  J.F.,  Thermal  conductivity  of  earth  materials  at 
high  pressures,  Ph.D.  thesis,  Massachusetts  Institute 
of  Technology,  1971. 


120 


Sclater,  J.G.,  and  Jean  Francheteau,  The  implications  of 

terrestrial  heat  flow  observation  on  current  tectonic 
and  geochemical  models  of  the  crust  and  upper  mantle 
of  the  earth,  Geophys.  J.  R.  astr.  Soc.,  20,  509, 

1970. 

Sleep,  N.H.,  Sensitivity  of  heat  flow  and  gravity  to  the 
mechanism  of  sea-floor  spreading,  J .  Geophys .  Res . , 

74,  542,  1969. 

Taiwan!,  M. ,  Gravity,  The  Sea,  Vol.  4,  edited  by  A.E. 
Maxwell,  Interscience,  New  York,  251,  1970. 

Taiwan!,  M. ,  X.  Le  Pichon,  and  M.  Ewing,  Crustal  structure 
of  the  mid-ocean  ridges,  computed  model  from  gravity 
and  seismic  refraction  data,  J.  Geophys.  Res.,  70, 

341,  1965. 

Taiwan!,  M. ,  C.C.  Windisch,  and  M.G.  Langseth,  Jr., 

Reykjanes  ridge  crest:  a  detailed  geophysical  study, 

J.  Geophys.  Res. ,  76 ,  473,  1971. 

Torrance,  K.E.,  and  D.L.  Turcotte,  Thermal  convection  with 
large  viscosity  variations.  Journal  Flu.  Mech. ,  47 , 
1971a. 

Torrance,  K.E.  and  D.L.  Turcotte,  Structure  of  convection 

cells  in  the  mantle,  J.  Geophys.  Res.,  76 ,  3154,  1971b. 

Varga,  R.S.,  Matrix  Iterative  Analysis,  Prentice-Hall,  New 


York,  1962. 


-  121  - 


Figure  1 


Figure  2. 


Figure  3. 


Figure  4. 


Figure  5. 


FIGURE  CAPTIONS 

Diagram  showing  grid  locations  to  which 
subscripts  refer  in  equations  (15)  and  (16), 
relative  to  grid  point  i,j  (point  0). 

Grid  locations  to  which  subscripts  refer 
in  equations  (17),  (18),  and  (19). 

Left:  Normal  geotherm  used  in  the  calculations. 
Center:  Viscosity  model  corresponding  to  the 
normal  geotherm.  Right:  Velocity  profile  of 
plate  motion  and  counterflow  below  that  mini¬ 
mizes  shear  dissipation  for  horizontal  flow. 

This  velocity  profile  was  used  at  the  artificial 
lateral  boundary  in  the  calculations. 

Streamlines  and  isotherms  for  plate  velocity  of 
1.2  cm/yx  in  steady-state. 

Heat  flow  as  a  function  of  age.  Average  heat 
flow  measurements  for  the  North  Pacific  [Sclater 
and  Francheteau.  1970 ]  are  compared  with 
calculated  heat  flow  for  Reyk janes  60°N,  North 
Atlantic  32°,  South  Atlantic  28°,  and  the  East 
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Figure  6. 


Figure  7. 


Figure  8. 


Pacific  43°S.  Error  bars  denote  standard 
deviation  of  obssrvsd  heat  flow  data. 

North  Atlantic  32°,  topography  and  gravity. 
Top:  Observed  Bouguer  anomaly  [Taiwan!  and 
Le  Pichon ,  1965]  is  compared  to  present 
calculation.  Bottom:  Observed  topography 
[Taiwan! ,  1970]  is  compared  to  calculated 
isostatic  topography  (dashed  curve)  and 
calculated  topography  resulting  from  viscous 
stresses  (solid  curve) .  Sea  floor  spreading 
rates  wera  averaged  from  Pitman  and  Talwani 
[1971]. 

South  Atlantic  28° ,  topography.  Observed 
topography  [Talwani.  1970]  is  compared  to 
calculated  isostatic  topography  (dashed  curve) 
and  calculated  topography  resulting  from 
viscous  stresses  (solid  curve) .  Sea  floor 
spreading  rates  were  averaged  from  Maxwell 
[1969],  Dickson  et  al. [19681  and  Heirtzler  et 
al.  [1968]. 

Reykjanes  60#N,  topography  and  gravity.  Top: 
Observed  Bouguer  anomaly  [Talwani  et  al. ,  1971] 
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is  compared  to  present  calculation.  Bottom: 
Observed  topography  is  compared  to  calculated 
isostatic  topography  (dashed  curve)  and  calculated 
topography  resulting  from  viscous  stresses 
(solid  curve) .  Sea  floor  spreading  rates  were 
taken  from  Pitman  and  Talwani  [1971] . 

Figure  9.  South  Pacific  43°,  topography.  Observed 
topography  [Talwani ,  1970]  is  compared  to 
calculated  isostatic  topography  (dashed  curve) 
and  calculated  topography  resulting  from 
viscous  stresses  (solid  curve).  Sea  floor 
spreading  rates  were  obtained  from  Pitman 
et_al.  [1968]  and  Heirtzler  et  al.  [1968]. 

Figure  10.  Horizontal  stress  deviator  integrated  through 
the  thickness  of  the  plate  and  divided  by 
60  km,  /  s^  dy/60,  is  plotted  as  a  function 
of  distance  from  the  ridge  axis. 
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Derr,  J.  s..  Discrimination  of  earthquakes  and  explosions  by 
the  Rayleigh-wave  spectral  ratio.  Bull.  Seism.  Soc.  Am. ■ 
60,  1653-1668,  1970. 

Gilbert,  F.  and  D.  V.  Helmberger,  Generalized  ray  theory  for 
a  layered  sphere,  J.  Oeophys.  Res..  in  press,  1971. 

Helmberger,  D.  and  R.  A.  Wiggins,  Upper  mantle  structure  of 
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midwestern  United  States ,  J.  Geophys.  Res./  76 ,  3229- 
3245,  1971. 

Hirasawa,  T.,  Radiation  patterns  of  S  waves  from  underground 

nuclear  explosions,  J.  Geophys.  Res.,  76,  6440-6454,  1971. 

Hirasawa,  T.  and  M.  J.  Berry,  Reflected  and  head  waves  from  a 
linear  transition  layer  in  a  fluid  medium.  Bull.  Seism. 

Soc .  Am. ,  61,  1-26,  1971. 

Husebye,  E.  and  R.  Madariaga,  The  origin  of  precursors  to  core 
waves.  Bull.  Seism.  Soc.  Am.,  60,  939-952,  1970. 

Mack,  H. ,  Twenty-  to  eighty-second  period  characteristics  of 
nuclear  explosions  recorded  at  LASA,  in  Copies  of  Papers 
Presented  at  Woods  Hole  Conference  on  Seismic  Discrimina¬ 
tion  ,  1 ,  ARPA  Working  Paper,  July  20-23,  1970. 

Mack,  H.,  Multipa thing  of  Rayleigh  waves  generated  by  Milrow, 
in  Copies  of  Papers  Presented  at  Woods  Hole  Conference  on 
on  Seismic  Discrimination,  2,  ARPA  Working  Paper,  July  20- 
23,  1970. 

Toksdz,  M.  N.,  J.  W.  Minear,  and  B.  R.  Julian,  Temperature  field 
and  geophysical  effects  of  a  downgoing  slab,  J .  Geophys . 
Res.,  76,  1113-1138,  1971. 
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VIII. 

RESEARCH  ASSOCIATES  SUPPORTED  BY  CONTRACT  NO.  AF4 9/638-1763 

Below  is  a  brief  resume  for  each  of  the  twelve  Research 
Associates,  in  alphabetical  order,  supported  by  the  present 
contract  during  the  years  1966  to  1971.  Included  are  their 
periods  of  tenure,  their  research  interests  and  the  publica¬ 
tions  that  resulted  from  work  done  at  M.  I.  T.  Their  studies 
and  contributions  in  seismology  cover  a  broad  spectrum  of 
problems.  They  all  had  some  exposure  to  Large  Aperture 
Seismic  Arrays,  theoretical  techniques  and  seismic  data.  A 
number  of  these  scientists  are  fully  involved  with  seismic 
arrays  and  research  on  seismic  discrimination. 


Dudley  J.  Andrews 

Dr.  Andrews  has  been  a  Research  Associate  since  1970. 

He  has  been  working  on  numerical  calculations  of  wave  propa¬ 
gation,  crustal  movement,  and  mantle  flow.  Of  particular 
interest  is  his  work  on  the  shape  of  the  pressure  pulse  from 
an  underground  nuclear  explosion  in  the  nearly  elastic  range 
of  distances. 

Publications : 

Andrews,  D.  J.,  A  numerical  method  for  creep  deformation  of 
solids  (abstract) ,  EOS  Trans.  Amer.  Geophys.  Un.,  52, 
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347,  1971. 

Andrews,  D.  J.  and  S.  Shlien,  Calculated  attenuation  of  an 
underground  explosion  wave  in  the  nearly  elastic  range 
(abstract),  Earthquake  Notes,  42,  9,  1971. 

Andrews,  D.  J. ,  Numerical  simulation  of  sea-floor  spreading, 
submitted  to  J.  Geophys.  Res.,  1971. 


William  H.  Bakun 

Dr.  Bakun  was  a  Research  Associate  during  1970-1971. 

He  studied  the  crustal  structure  beneath  LASA  using  the 
Haskell -Thomson  theory  and  he  determined  the  focal  depths  of 
earthquakes  from  near-source  observations  of  various  crustal 
P  phases . 

Publications  include : 

Bakun,  W.  H.,  Focal  depths  of  the  1966  Parkfield,  California, 
earthquakes,  J.  Geophys.  Res.,  in  press,  1971. 

Bakun,  W.  H.,  Crustal  structure  beneath  LASA  from  long-period 
P-wave  spectra,  submitted  to  J.  Geophys.  Res.,  1971. 

David  M.  Boore 

Dr.  Boore  was  a  Research  Associate  during  1969-1970. 

His  interests  included  numerical  and  wave  theoretical  studies 
of  body  and  surface  wave  propagation  in  heterogeneous  media 
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media  and  theoretical  models  of  the  earthquake  source. 
Publications : 

Boore ,  D.,M.,  K.  Aki,  and  T.  Todd,  A  two-dimensional  moving 
dislocation  model  for  a  strike-slip  fault,  Bull.  Seism. 
Soc.  Am.,  61,  177-194,  1971. 

Boore,  D.  M,  K.  L.  Larner,  and  K.  Aki,  Comparison  of  two  inde¬ 
pendent  methods  for  the  solution  of  wave  scattering  pro¬ 
blems:  response  of  a  sedimentary  basin  to  vertically 
incident  SH  waves,  J.  Geophys.  Res.,  76,  558-569,  1971. 


John  S.  Derr 

Dr.  Derr  was  a  Research  Associate  during  1968-1970.  He 
analyzed  data  from  M.I.T.'s  new  mercury  tiltmeter  and  studied 
the  amplitude  spectra  of  Rayleigh  waves  as  a  discriminant  be¬ 
tween  earthquakes  and  explosions. 

Publications : 

Derr,  J.  S.,  Discrimination  of  earthquakes  and  explosions  by 
the  Rayleigh  wave  spectral  ratio.  Bull.  Seism.  Soc.  Am., 
60,  1653-1668,  1970. 


Donald  V.  Helmberger 

Dr.  Helmberger  was  a  Research  Associate  during  1967-1969. 
He  worked  on  generalized  ray  theory  and  its  use  in  the  compu- 
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tation  of  theoretical  seismograms.  By  comparing  theoretical 
with  observed  records,  he  studied  the  compressional  and  shear 
wave  velocities  of  the  crust  and  upper  mantle  and  the  fine 
structure  of  the  crust-mantle  transition. 

Publications  include : 

Helmberger,  D.  V.  and  G.  B.  Morris,  A  travel  time  and  amplitude 
interpretation  of  a  marine  refraction  profile:  primary 
waves'  Geophys.  Res. ,  74.  483-494,  1969. 

Helmberger,  D.  V.  and  G.  B.  Morris,  A  travel  time  and  amplitude 
interpretation  of  a  marine  refraction  profile:  transformed 
shear  waves.  Bull.  Seism.  Soc.  Am.,  60.  593-600,  1970. 
Helmberger,  D.  V.  and  R.  A.  Wiggins,  Upper  mantle  structure 
of  midwestern  United  States,  J.  Geophvs.  Res. .  3229- 

3245,  1971. 

Gilbert,  F.  and  D.  v.  Helmberger,  Generalized  ray  theory  for 
a  layered  sphere,  in  press  J.  Geophvs.  Res..  1971. 


Tomowo  Hirasawa 

Dr.  Hirasawa  was  a  Research  Associate  during  1968-1969. 

His  work  included  wave  theoretical  studies  of  compressional 
waves  and  an  investigation  of  tectonic  strain  release  by  nuclear 
explosions  using  shear  wave  radiation  patterns. 

Publications  include: 

Hirasawa,  T.  and  M.  J.  Berry,  Reflected  and  head  waves  from  a 
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linear  transition  layer  in  a  fluid  medium.  Bull.  Seism.  Soc. 

Am. ,  61,  1-26,  1971. 

Hirasawa,  T. ,  Radiation  patterns  of  S  waves  from  underground 
nuclear  explosions,  J.  Geophys.  Res.,  76,  6440-6454,  1971. 


Eystein  S.  Husebye 

Dr.  Husebye  was  a  Research  Associate  during  1960-1968. 

He  was  interested  in  the  seismic  array  as  a  tool  for  studying 

the  earth.  He  developed  methods  for  measuring  dT/dA  and 
2  2 

d  T/dA  using  LASA  and  applied  the  techniques  to  study  core 
phases.  Husebye  is  currently  on  the  staff  of  NORSAR. 
Publications  include: 

Husebye,  E.  S.,  Direct  measurement  of  dT/dA,  Bull.  Seism. 
Soc.  Am. ,  59,  717-727,  1969. 

Husebye,  E.  and  R.  Madariaga,  The  origin  of  precursors  to  core 
waves,  Bull.  Seism.  Soc.  Am.,  60,  939-952,  1970. 

Toksdz,  M.  N.,  E.  S.  Husebye,  and  R.  A.  Wiggins,  Structure  and 
properties  of  the  earth's  core,  in  preparation,  1971. 


Bruce  R.  Julian 

Dr.  Julian  was  a  Research  Associate  during  1969-1970. 

His  work  focused  on  ray  theory  for  heterogeneous  media  with 
applications  to  lateral  variations  of  seismic  velocity  in  the 
mantle.  Particular  emphasis  was  placed  on  understanding  the 
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effect  on  body-wave  travel -times  and  amplitude  radiation  pat¬ 
terns  of  downgoing  plates.  Dr.  Julian  is  currently  a  member 
of  the  Seismic  Discrimination  Group  at  M.  I.  T.  Lincoln  Labora¬ 
tory. 

Publications  include : 

Julian,  B.  R. ,  Regional  variations  in  upper  mantle  structure  in 
North  America  (abstract) ,  EOS  Trans.  Amsr.  Geophys.  Un., 
51,  359,  1970. 

Toksdz,  M.  N. ,  J.  w.  Minear,  and  B.  R.  Julian,  Temperature 
field  and  geophysical  effects  of  a  downgoing  slab,  J. 
Geophys.  Res.,  76,  1113-1138,  1971. 


Harry  Mac); 

Dr.  Mack  was  a  Research  Associate  during  1969-1970.  He 
worked  on  array  processing  of  seismic  and  microbarograph  data. 

In  particular,  using  freqeuncy-wave  number  analysis,  he  studied 
the  spectral  characteristics  of  Rayleigh  waves  from  several 
nuclear  explosions  and  earthquakes  and  the  problem  of  surface- 
wave  multipathing.  Dr.  Mack  is  currently  with  the  Array  Analysis 
Center  at  Teledyne. 

Publications  include: 

Mack,  H.,  Twenty-  to  eighty-second  period  characteristics  of 
nuclear  explosions  recorded  at  LASA,  in  Copies  of  Papers 
Presented  at  Woods  Hole  Conference  on  Seismic  Discrimina- 
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tj-on'  1  (ARPA) ,  July  20-23,  1970. 

Mack,  H.,  Multipa thing  of  Rayleigh  waves  generated  by  Milrow, 
in  Coeies  of  Papers  Presented  at  Woods  Hole 
on.  Seismic  Discriminat j on r  7  (ARPA),  July  20-23,  1970. 

John  w.  Minear 

Dr.  Minear  was  a  Research  Associate  during  1968-1969. 

He  constructed  numerical  models  for  mantle  convection  and 
lithosphere  subduction.  Particularly  important  work  was  his 
temperature  models  for  downgoing  slabs  of  lithosphere  and  the 

associated  effects  on  such  geophysical  properties  as  seismic 
wave  propagation. 

Publications  include: 

Minear,  J.  w.  and  M.  N.  ToksSz,  Thermal  regime  of  a  downgoing 
slab  and  new  global  tectonics,  J.  Geophvs.  Re...  75, 
1397-1419,  1970. 

Minear,  J.  W.  and  M.  M.  Toksdz,  Thermal  regime  of  a  downgoing 
Slab'  Tectonophysics,  10.  367-390,  1970. 

Toksdz,  M.  N.,  J.  w.  Minear,  and  B.  R.  Julian,  Temperature 
field  and  geophysical  effects  of  a  downgoing  slab,  J^ 
Geophys.  Res.,  76,  1113-1138,  1971. 


Masanori  Saito 


Dr.  saito  was  a  Research  Associate  during  1966-1968.  His 
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work  included  a  technique  for  synthesizing  dilatational  and 
rotational  seismograms  from  an  array  of  three-component  sta¬ 
tions  and  the  computation  of  variational  parameters  and  exci¬ 
tation  functions  for  surface  waves  and  free  oscillations  in 
layered  media.  His  computer  programs  from  his  latter  studies 
have  been  used  extensively  by  students  and  staff  at  M.I.T.  and 
elsewhere  to  study  earth  structure  and  source  properties  from 
surface  wave  dispersion  and  amplitude  data. 

Publications  include: 

Saito,  M. ,  Excitation  of  free  oscillations  and  surface  waves 
by  a  point  source  in  a  vertically  heterogeneous  earth, 

J.  Geophys.  Res.,  72,  3689-3699,  1967. 

Saito,  M.,  Synthesis  of  rotational  and  dilatational  seismograms, 
J.  Phys.  Earth,  16,  53-61,  1968. 

Saito,  M. ,  Partial  derivatives  of  the  phase  velocity  of  surface 

waves  with  respect  to  anisotropy  factors,  in  preparation, 
1971. 


Roger  Turpening 

Dr.  Turpening  was  a  Research  Associate  during  1966-1967. 
He  studied  particle  motion-mode  filter  and  designed  a  digital 
filter  to  enhance  rectilinear  motion  in  the  presence  of  retro¬ 
grade  elliptical  motion. 


